Análisis estadístico de datos de Microbioma con R
Capítulo 7: Análisis exploratorio de datos de microbiomas
Equipo 4
Andrés Arredondo Cruz (andresabstract@gmail.com)
Adriana Haydé Contreras Peruyero (haydeeperuyero@gmail.com)
David Alberto García Estrada (david.garcia.e@cinvestav.mx)
2022-10-04
Datos de Ratones y Humanos
Conjuto de datos de ratón \(Vdr^{-/-}\)
Los datos del microbioma intestinal murino (Jin et al.2015) se obtuvieron de muestras de heces fecales y cecales. Aquí, se utilizan las muestras fecales.
Conjuto de datos de fumadores de cigarro
El segundo conjunto de datos es de Charlson et al. (2010) y Chen (2012), que incluye estudios sobre el efecto del tabaquismo en el microbioma del tracto respiratorio superior. El conjunto de datos original contiene muestras de microbiomas de garganta y nariz, y de ambos lados del cuerpo. El conjunto de datos utilizado en este capítulo proviene del microbioma de la garganta del lado izquierdo del cuerpo. Contiene 60 sujetos (32 no fumadores y 28 fumadores). El conjunto de datos incluye tres datos: recuento de abundancia, árbol y metadatos. Es adecuado para ilustrar el diagrama de árbol y el análisis de ordenación restringida.
Análisis exploratorio como resumen gráfico
Podemos dividir los métodos de estudio de la composición de la comunidad de microbiomas en dos componentes principales: - Análisis de diversidades taxonómicas - Análisis multivariante de la composición del microbioma
El análisis multivariado incluye varias técnicas multivariadas, como el agrupamiento y la ordenación (sin restricciones y con restricciones) y las diferencias de prueba de hipótesis entre los grupos.
Usaremos el paquete de “Phyloseq” que es una herramienta bastante importante, almacena, analiza y muestra gráficamente datos complejos filogenéticos. Los datos de entrada usados en este paquete pueden ser OTUs o datos de abundancia de cuentas. Este paquete usa sistemas gráficos y avanzados (ggplot2) para facilitar gráficas de calidad de datos.
En este capítulo se exploran diferentes gráficos usuales: riqueza, barras de abundancia, mapas de calor, redes y árbol filogenético.
Gráficos de riqueza
Las diversidades alfa estimadas pueden resumirse mediante un gráfico
utilizando la función plot_richness() del paquete
phyloseq. Aunque su nombre sugiere trazar “riqueza”, que
normalmente se refiere a trazar el número total de especies/taxa/OTUs en
una muestra o entorno, en realidad la función no sólo traza la riqueza,
sino que también genera cifras de diversidades observadas y otras
estimadas.
En primer lugar, debemos cargar los paquetes phyloseq y
ggplot2, y el conjunto de datos de ratones Vdr\(-/-\).
#Para instalar el paquete, usamos Bioconductor
#if (!require("BiocManager", quietly = TRUE))
# install.packages("BiocManager")
#BiocManager::install("phyloseq")# IMPORTANTE
# Ajustamos path de trabajo según tu PC
# Mostramos directorio actual
getwd()## [1] "D:/Users/hayde/Documents/R_sites/Equipo4/Chapter7"# Si es necesario, cambiar la ruta donde están los archivos almacenados.
# "." significa el directorio actual.
workingDir <- "."
setwd(workingDir)# cargamos librerias
library("phyloseq")
library("ggplot2")
library("igraph")##
## Attaching package: 'igraph'## The following objects are masked from 'package:stats':
##
## decompose, spectrum## The following object is masked from 'package:base':
##
## unionlibrary("vegan")## Loading required package: permute##
## Attaching package: 'permute'## The following object is masked from 'package:igraph':
##
## permute## Loading required package: lattice## This is vegan 2.6-2##
## Attaching package: 'vegan'## The following object is masked from 'package:igraph':
##
## diversitylibrary("GUniFrac")
library("pbkrtest")## Loading required package: lme4## Loading required package: Matrixlibrary("BiodiversityR")## Loading required package: tcltk## BiodiversityR 2.14-3: Use command BiodiversityRGUI() to launch the Graphical User Interface;
## to see changes use BiodiversityRGUI(changeLog=TRUE, backward.compatibility.messages=TRUE)library("kableExtra")## Registered S3 method overwritten by 'httr':
## method from
## print.response rmutilEn la carpeta data se encuentra la base de datos a usar
VdrFecalGenusCounts.csv.
# Mandamos llamar los datos originales
# Fila es el taxón y columna la condición
abund_table=read.csv(paste0(workingDir,"/data/VdrFecalGenusCounts.csv"), row.names=1,check.names=FALSE)
# Generamos la transpuesta, fila la condición, columna el taxón
abund_table<-t(abund_table)El paso crítico cuando se utiliza el paquete phyloseq es construir un
objeto phyloseq-class. El siguiente código R construye un
objeto de clase phyloseq llamado physeq utilizando el comando
phyloseq(). El objeto de clase phyloseq se construye a
partir de sus componentes datos: - Tabla OTU. - Datos muestra. - Tabla
de taxonomía. - Árbol filogenético.
Como objeto a nivel de experimento, deben proporcionarse dos o más objetos de datos componentes. El orden de los argumentos no importa.
# Hacemos una tabla que incluye el nombre de las filas y además incluye como dato,
#estos nombres separados por "_"
meta_table <- data.frame(row.names=rownames(abund_table),
t(as.data.frame(strsplit(rownames(abund_table),"_"))))
# Hacemos factor la columna "X2" y los que sean 11, 12, 13, 14 y 15 serán "Vdr-/-",
#el resto será "WT"
meta_table$Group <- with(meta_table,ifelse(as.factor(X2) %in%
c(11,12,13,14,15), c("Vdr-/-"), c("WT")))Convertimos los datos al formato phyloseq.
# OTUs
OTU <- otu_table(as.matrix(abund_table), taxa_are_rows = FALSE)
# Sample
SAM <- sample_data(meta_table)
# Unimos en un objeto phyloseq
physeq <- merge_phyloseq(phyloseq(OTU),SAM)
physeq## phyloseq-class experiment-level object
## otu_table() OTU Table: [ 248 taxa and 8 samples ]
## sample_data() Sample Data: [ 8 samples by 4 sample variables ]Una vez que tenemos nuestro objeto phyloseq, podemos
usar la función plot_richness() para construir graficar las
diversidades alpha observadas y estimadas.
plot_richness(physeq, x = "Group ", color = "Group ")
Gráficos de diversidad alpha con Vdr y grupos WT en muestras de heces.
Esta función también nos permite seleccionar solo algunas diversidades. El siguiente es un ejemplo usando solo dos diversidades, la de Chao1 y Shannon.
plot_richness(physeq, measures = c("Chao1", "Shannon"),x = "Group ", color = "Group ")
Gráficos de diversidad alpha seleccionando las diversidades de Chao y Shannon.
Se requiere el dato de entrada “physeq”, que es de la clase
“phyloseq”, o alternativamente, un una tabla OTUs. El argumento opcional
“x” es una variable que se asigna al eje horizontal; x puede ser una
cadena de caracteres o un vector. El valor por defecto es “samples”, que
asignará el nombre de cada muestra a una posición horizontal distinta en
el gráfico. En este caso, x = "Group" asignará la
pertenencia a un grupo al eje x. El argumento color también es opcional.
Especificará la variable de muestra que se asignará a diferentes
colores. Al igual que el argumento x, puede ser una cadena de caracteres
o un vector. En la estimación de la diversidad alfa, el ruido puede
recortarse; sin embargo, dado que muchas estimaciones de riqueza e
incluso la riqueza “observada” dependen en gran medida del número de
individuos. Por lo tanto, si desea obtener resultados significativos,
debe utilizar conjuntos de datos sin recortar (McMurdie y Holmes
2013).
Gráfico de barras de abundancia
El diagrama de barras por defecto, sin ningún parámetro, trazará con cada muestra individualmente en el eje x, y los valores de abundancia en el eje y. Los valores de abundancia para cada OTU/muestra se apilan en el orden de mayor a menor, separados por una fina línea horizontal.
theme_set(theme_bw())
# Gráfico de barras
plot_bar(physeq)
Gráfico de barras de abundancia por defecto.
Otros parámetros que se pueden proporcionar son fill y
facet_grid. El parámetro fill es una cadena de
caracteres que especifica cual variable muestral puede ser usada para
mapearla a los colores con los cuales se llenan las barras.
plot_bar(physeq, fill="Group")
Gráfico de barras de abundancia indicando cual variable es mapeada al parámetro fill.
El parámetro facet_grid es una fórmula que especifica
las facetas a mostrar en el gráfico. En el siguiente código, primero
encontramos los 5 taxones con más abundancia y después los graficamos
añadiendo también el parámetro de facet_grid.
# Nos da el nombre de los taxones top 5
TopNGenus <- names(sort(taxa_sums(physeq), TRUE)[1:5])
# Generamos un objeto phyloseq con solo los 5 taxa que escogimos
Top5Genus <- prune_taxa(TopNGenus, physeq)
# Graficamos
plot_bar(Top5Genus, fill="Group", facet_grid=~Group)
Gráfico de barras de abundancia de los 5 taxones con mayor abundacia.
Otra opción que podríamos especificar dentro del parámetro
fill es Genus pero estos datos no cuentan con
ese parámetro.
# plot_bar(Top5Genus, fill="Genus", facet_grid=~Group)Abundancias relativas
Calculamos las abundancias relativas con la función \(x/sum(x)*100\). Después convertimos nuestro
objeto top_5_relativo a un data frame para poder usar como
factores los OTU y muestras y de esta forma poder llenar por OTU nuestro
gráfico de barras.
#Abundancias relativas
top_5_relativo <- transform_sample_counts(Top5Genus, function(x) x/sum(x)*100)
#Grafico de barras con abundancias relativas
plot_bar(top_5_relativo, fill="Group", facet_grid = ~Group)
Gráfico de barras de abundancia de los 5 taxones con mayor abundacia relativa.
#Convertimos a dataframe
top_5_relativo_df <- psmelt(top_5_relativo)
#Damos formato como caracter y factor a los OTU
top_5_relativo_df$OTU <- as.character(top_5_relativo_df$OTU)
top_5_relativo_df$OTU <- as.factor(top_5_relativo_df$OTU)
#Convertimos a factor las muestras, no fue necesario determinar los niveles, ya los
#detecta como tal
top_5_relativo_df$Sample <- as.factor(top_5_relativo_df$Sample ) #, levels =c("5_15_drySt-28F", "1_11_drySt-28F", "2_12_drySt-28F", "3_13_drySt-28F", "4_14_drySt-28F", "7_22_drySt-28F", "8_23_drySt-28F", "9_24_drySt-28F"))
#Damos formato con ggplot
library(RColorBrewer)
#definimos paleta de colores por niveles de OTU
OTU_colors<- brewer.pal(length(levels(top_5_relativo_df$OTU)),"Dark2")
#especificamos los aes de ggplot
global_rel_plot<-ggplot(top_5_relativo_df, aes(x=Sample, y=Abundance, fill=OTU)) +
geom_bar(aes(), stat="identity", position="stack") +
scale_fill_manual(values=OTU_colors) +
labs(x= "Sample", y = "Relative abundance") +
facet_grid(~Group, scales = "free_x")+
ggtitle("Relative abundance by Groups")+
theme(axis.title = element_text(size = 14), axis.text.x = element_text(size = 5), plot.title = element_text(size = 20))
global_rel_plot
Gráfico de barras de abundancia de los 5 taxones con mayor abundacia relativa por OTU.
Gráfico de mapa de calor
El ordenamiento basado en la ordenación es mucho mejor que el cluster
jerárquico para presentar los datos del microbioma. Aquí, nos centramos
en ilustrar la función plot_heatmap() basada en los métodos
de ordenación NMDS y PCA. A continuación se presenta un uso de la
función plot_heatmap().
plot_heatmap(physeq_object, method = “NMDS”, distance = “bray”, sample.label = NULL, taxa.label = NULL, low = “#000033”, high = “#66CCFF”, na.value = “black”).
Se requiere el argumento de datos de entrada “physeq”. Es un objeto
de la clase phyloseq (otu_table). Tanto el método como los aumentos de
distancia son opcionales. El método de ordenación es el que se utilizará
para organizar el mapa de calor. El método de distancia ecológica es una
cadena de caracteres para usar en la ordenación. Tanto “sample.label”
como “taxa.label” son cadenas de caracteres y opcionales para usar para
etiquetar el eje de la muestra (horizontal) y reetiquetar el eje de
taxones/especies/OTU (vertical), respectivamente. Los aumentos bajos y
altos son cadenas de caracteres y opcionales. Se utilizan para elegir
las opciones de color compatibles con R. R entiende más de 600 colores.
Puede escribir colors() en R para comprobar los
nombres.
En los mapas de calor, los valores \(0\) son tratados como valores faltantes (NA) y se mapean al color negro. Por lo general, los valores más pequeños se representan por default con el color azul obscuro; mientras que los valores altos con los colores más claros. Para ejemplificar esto, vamos a usar primero solo los cinco taxones con mayor abundancia.
# Nos da el nombre de los taxones top 5
TopNGenus <- names(sort(taxa_sums(physeq), TRUE)[1:5])
# Generamos un objeto phyloseq con solo los 5 taxa que escogimos
Top5Genus <- prune_taxa(TopNGenus, physeq)
# Graficamos
plot_heatmap(Top5Genus) 
Mapa de calor de los 5 taxones con mayor abundancia.
plot_heatmap(top_5_relativo)
Mapa de calor de los 5 taxones con mayor abundancia relativa.
Vamos a especificar otros parámetros como la distancia y métodos.
# Método NMDS y distancia de Bray-Curtis
p <- plot_heatmap(Top5Genus, "NMDS", "bray")
p
Mapa de calor de los 5 taxones con mayor abundancia usando el métrodo NMDS y la distancia Bray-Curtis.
Podemos especificar también los rangos de colores.
plot_heatmap(Top5Genus, "NMDS", "bray", low="#000033", high="#CCFF66") 
Mapa de calor especificando otras escalas de colores.
plot_heatmap(Top5Genus, "NMDS", "bray", low="#000033", high="#FF3300")
Mapa de calor especificando otras escalas de colores.
plot_heatmap(Top5Genus, "NMDS", "bray", low="#000033", high="#66CCFF")Mapa de calor especificando otras escalas de colores.
plot_heatmap(Top5Genus, "NMDS", "bray", low="#66CCFF", high="#000033", na.value="white")
Mapa de calor especificando otras escalas de colores.
El siguiente mapa de calor es un ejemplo usando PCoA.
# Método de PCoA y distancia de Bray-Curtis
plot_heatmap(Top5Genus, "PCoA", "bray")
Mapa de calor usando PCoA.
Gráfico de redes
Hay dos funciones en el paquete phyloseq para trazar la red del
microbioma usando “ggplot2”: plot_network() y
plot_net(). Si se utiliza la función
plot_network(), su primer argumento debe ser un objeto
igraph, y la red en sí debe representarse utilizando el
paquete igraph. El objeto red (argumento g) puede ser creado usando la
función make_network() a través del paquete phyloseq.
La función plot_net() es una revisión del rendimiento y
de la interfaz de plot_network(), su primer/principal
argumento es una instancia de la clase phyloseq. Los ejemplos de uso de
plot_network() y plot_net() se dan a
continuación:
plot_network(g, physeq=NULL, type=“samples”, color=“Group”, shape=“Group”)- g es un objeto de clase igraph necesario creado por la función
make_network(), o directamente por el paquete igraph. - El argumento opcional physeq es un objeto de clase phyloseq en el que se basa g.
- La opción
Typeindica si el grafo representado en el argumento primario, g, es de muestras o de taxones/OTUs. Por defecto es “samples”. - Las opciones color y shape son opcionales. Se utilizan para el mapeo de color y el mapeo de forma de los puntos.
- g es un objeto de clase igraph necesario creado por la función
plot_net(physeq, distance = “bray”, type = “samples”, maxdist = 0.7, color = NULL, shape = NULL).- El argumento requerido physeq, es el objeto de clase phyloseq que se quiere representar como una red.
- La opción de distancia es un método de distancia o una clase-distancia ya calculada. Por defecto es bray.
- La opción
maxdistsignifica el valor máximo de distancia entre dos vértices para conectar con una arista en el grafo. El valor por defecto es \(0,7\).
Los siguientes códigos de R crean un grafo basado en igraph, basado en el método de distancia por defecto, Jaccard y una distancia máxima entre nodos conectados de \(0,8\). El “Group” se utiliza para los mapeos de color y forma para visualizar la estructura de las muestras de ratones Vdr\(-/-\) y WT.
set.seed(123)
library("igraph")
# Hacemos una gráfica a partir de el objeto phyloseq
ig <- make_network(physeq, max.dist=0.8)
# Graficamos
plot_network(ig, physeq, color="Group", shape="Group")
Gráfico usando plot network().
El gráfico anterior muestra una estructura interesante, con dos subgráficos que comprenden las muestras de los ratones Vdr-/- y WT respectivamente. Además, parece parece haber una correlación entre las muestras.
En comparación con la función plot_network(), la nueva
función plot_net() no requiere una llamada separada a la
función make_network(), o un objeto igraph separado. Los
códigos siguientes códigos crean una red basada en una distancia máxima
entre los nodos conectados de \(0,5\).
Como estamos interesados en cómo la estructura de las muestras de
ratones Vdr-/- y WT muestras, utilizamos el “Group”
para el mapeo de color y el mapeo de forma de los puntos.
# Graficamos sin necesidad de hacer el objeto anterior
plot_net(physeq, maxdist = 0.5, color = "Group", shape="Group")
Gráfico de red creado con plot net().
Gráfica de árbol filogenético
La función plot_tree() del paquete phyloseq pretende
facilitar la investigación gráfica del árbol filogenético, así como de
los datos de la muestra. Para la secuenciación filogenética de muestras
con gran riqueza, el árbol generado por esta función será difícil de
leer e interpretar. Una “regla general” propuesta por los autores del
paquete phyloseq es utilizar subconjuntos de datos con no más de 200 OTU
por parcela. A continuación se muestra un uso de esta función:
plot_tree(physeq, method = “sampledodge”, color = NULL, shape = NULL, ladderize = FALSE)
Los datos de entrada “physeq” son necesarios. Los datos deben ser de la clase phyloseq y contener como mínimo un componente de árbol filogenético. Son estos datos los que se quieren trazar y anotar un árbol filogenético. La opción “method” es el nombre del método de anotación a utilizar. El método por defecto “sampledodge” el cual dibujará puntos junto a las hojas si se observaron individuos de ese taxón, y un punto separado para cada muestra. Los argumentos “color” y “forma” son opcionales. Proporcionan los nombres de las variables en physeq para mapear el color del punto, y mapear la forma del punto, respectivamente. La opción “ladderize” se utiliza para especificar si el árbol se ordena en forma de escalera o no, es decir reordenar los nodos de acuerdo con la profundidad de sus subárboles antes de trazarlos. Por defecto es FALSE, no se aplica la jerarquización. Cuando es TRUE o “right”, se utiliza la se utiliza la jerarquización. Cuando se establece como “izquierda”, se aplica la escalonación “izquierda”. El conjunto de datos de los fumadores se utiliza para ilustrar el trazado del árbol filogenético. Los datos son del paquete GUniFrac (Charlson et al. 2010; Chen 2012). Primero vamos a cargar este paquete y hagamos que los datos estén disponibles para su uso.
library("GUniFrac")
data(throat.otu.tab)
data(throat.tree)
data(throat.meta)Ahora, necesitamos construir nuestro objeto phyloseq.
# Convertimos los datos en un objeto phyloseq
# Tabla de OTUs
OTU <- otu_table(as.matrix(throat.otu.tab), taxa_are_rows = FALSE)
# Tabla de Sample
SAM <- sample_data(throat.meta)
# Árbol filogenético
TRE <- throat.tree
# Generación del objeto de phyloseq
physeq <- merge_phyloseq(phyloseq(OTU), SAM, TRE)Verificamos cuantos taxones/OTUs tienen nuestros datos.
# Número de taxones
ntaxa(physeq)## [1] 856Y seleccionamos solo los primeros \(50\).
# Nos quedamos con los primeros 50
physeq <- prune_taxa(taxa_names(physeq)[1:50], physeq) Ahora, graficamos usando como variable en el color el estado de fumador.
plot_tree(physeq, ladderize = "left", color = "SmokingStatus") 
Árbol filogenético usando como color la variable SmokingStatus.
plot_tree(physeq, ladderize = "left", color = "SmokingStatus", shape = "SmokingStatus") 
Árbol filogenético usando como color y forma la variable SmokingStatus.
Otro tipo de gráfico que se suele ocupar es al mapear los datos a coordenadas polares.
plot_tree(physeq, ladderize = "left", color = "SmokingStatus", shape = "SmokingStatus") +
coord_polar(theta = "y")
Árbol filogenético con coordenadas polares
Agrupamiento (Clústers)
Introducción a la agrupación, la distancia y la ordenación
La agrupación (o clasificación) y la ordenación son las dos clases principales de métodos multivariantes que suelen emplear los investigadores del microbioma y los ecologistas de comunidades. Hasta cierto punto, estos dos enfoques son complementarios. El objetivo de la agrupación es clasificar las muestras en clases (quizás jerárquicas) para reducir la complejidad (dimensionalidad) de los datos; puede reducir todas las muestras en una dimensión (eje x).
Sin embargo, los datos con dos o tres dimensiones son más interpretables que los con una dimensión porque la mayoría de los datos comunitarios son continuos. Con la ordenación los datos comunitarios pueden reducirse a dos o tres dimensiones. Por lo tanto, la ordenación suele ser deseada por los investigadores del microbioma y los ecologistas de comunidades. Se han desarrollado muchos métodos multivariantes basados en la ordenación en el estudio del microbioma y la ecología. Cuando se realiza la agrupación y la ordenación, es necesario que un método de distancia proporcione una medida de distancia.
Distancias y disimilitudes
Antes de entrar al tema de análisis de clústers es importante introducir los conceptos de distancias y disimilitudes o similitudes.
Una medida de disimilitud o distancia entre dos sujetos muy distintos será grande y por el contrario, una medida de similitud entre ellos será pequeña.
Una distancia es una función \(d: \mathcal{X}\times \mathcal{X}\to \mathbb{R}\) que cumple las siguientes propiedades (ver ): Dos de las distancias más comunes en estadística para variables aleatorias multivariadas cuantitativas son la distancia euclidiana y la distancia de Mahalanobis:Existen otras distancias, como la distancia de Gower que puede ser utilizada para estudiar variables aleatorias con entradas cuantitativas y/o cualitativas.
Informalmente, una similitud es una medida entre dos objetos de que tanto se parecen. Por lo cual entre más parecidos sean los dos objetos su medida de similitud será más grande. Usualmente toma valores positivos entre 0 y 1.
Formalmente, una similitud es una función \(s: \mathcal{X}\times \mathcal{X} \to \mathbb{R}\) que típicamente satisface las siguientes dos condiciones:Usualmente, la función de similitud satisface \(s(x,y)\geq 0\) pero existen ejemplos de similitudes que no lo satisfacen, por ejemplo: el coeficiente de correlación y los productos escalares. Sin embargo, si la función de similitud no es positiva y está acotada siempre es posible transformarla en una función de similitud positiva añadiendo una constante, es decir, podemos escribir nuestra nueva función de similitud \(\tilde{s}\) como \(\tilde{s}(x,y)=s(x,y)+c\), donde \(c\) es una constante positiva lo suficientemente grande (ver y ).
A diferencia de una distancia o una disimilitud, no existe un análogo para la desigualdad del triángulo para la similitud.
Informalmente, una disimilitud entre dos objetos es una medida del grado en el cual dos objetos son diferentes. Como es de esperarse, una disimilitud debe ser baja si se consideran pares de objetos que son similares.
De manera precisa, una disimilitud es una función \(d: \mathcal{X}\times \mathcal{X} \to \mathbb{R}\) que satisface las siguientes dos condiciones:Es importante mencionar que, en la literatura de análisis de datos como data mining y algunos textos de estadística (ver y ), a las disimilitudes también se les nombra distancias, lo cual es incorrecto. Lo que si ocurre es que toda distancia (o métrica) es una disimilitud, pero no toda disimilitud es una distancia.
En se menciona que las disimilitudes usuales de acuerdo al tipo de dato que se esté trabajando son las siguientes:
Tipo nominal: La disimilitud más sencilla es \[ d(x,y) = \left\{ \begin{matrix} 0 & \text{ si } x= y \\ 1 & \text{ si } x\neq y \end{matrix} \right. \]
Tipo ordinal: La disimilitud más sencilla es \[ d(x,y) = \frac{|x-y|}{n-1} \] donde suponemos que hay \(n\) datos, y a cada uno de ellos se les ha asignado un valor entero de \(0\) a \(n-1\).
Tipo intervalo: La disimilitud más sencilla es \[ d(x,y) = \| x- y\| \] donde \(\| \cdot \|\) denota la distancia euclidiana usual.
Distancias y disimilitudes entre datos
En el Cuadro siguiente, se muestra cómo se calculan distancias y disimilitudes entre individuos con datos numéricos, binarios, categóricos y mixtos.
tabla de distancias
Agrupamiento (Clústers)
Existen varias familias de métodos de agrupamiento disponibles en la literatura (Legendre y Legendre 2012): algoritmos secuenciales o algoritmos simultáneos, aglomerativos o divisivos, monotéticos o politéticos, métodos jerárquicos o métodos no jerárquicos, métodos probabilísticos frente a métodos no probabilísticos. Entre estas categorías, los métodos de agrupación jerárquica de Ward y de partición de k-means son los más comunes en los estudios sobre el microbioma.
Los métodos jerárquicos aglomerativos, dan origen a un gráfico llamado dendrogramas. Estos métodos son iterativos y en cada paso debe recalcularse la matriz de distancias, lo cual para bases de datos muy grandes suele consumir mucho tiempo de cómputo.
En este capítulo se ilustran varios métodos de clustering diferentes, incluyendo:
- Single linkage agglomerative clustering o método de la liga sencilla o del vecino más cercano.
- Complete linkage agglomerative clustering o método de la liga completa o del vecino más lejano.
- Average linkage agglomerative clustering o métdod de la liga promedio.
- Ward’s minimum variance clustering o método Ward.
El conjunto de datos de ratones Vdr\(-/-\) con muestras fecales con el que se ha estado trabajando es el que se utilizará para ilustrar el análisis de conglomerados. Aquí, la distancia Bray-Curtis se utilizará para ilustrar la clasificación de las muestras. También se aplican otras distancias.
Cargamos el paquete Vegan. Primero vamos a
normalizar la tabla de abundancia utilizando la función
decostand() y calcular las disimilitudes
Bray-Curtis entre todos los pares de muestras utilizando la
función vegdist() del paquete Vegan.
library("vegan")
abund_table_norm <- decostand(abund_table, "normalize")
bc_dist <- vegdist(abund_table_norm , method = "bray")Clustering aglomerativo de la liga sencilla
“Aglomerar” significa reunir en un clúster; la agrupación aglomerativa de la liga sencilla también se denomina clasificación de vecinos más cercanos. En cada paso, la agrupación combina dos muestras que contienen las distancias más cortas entre pares (o la mayor similitud), es decir, se toma la mínima de todas las posibles distancias entre los objetos que pertenecen a distintos conglomerados:
\[d_{AB}=\min_{i\in A, j\in B} (d_{ij}).\]
tabla de distancias
El resultado de la agrupación puede visualizarse como un dendrograma. El inconveniente del dendrograma resultante de un clustering de liga sencilla a menudo muestra el encadenamiento de las muestras. Los clústeres formados se ven forzados a unirse debido a que los elementos individuales están cerca unos de otros, mientras que muchos elementos de cada clúster pueden estar muy alejados entre sí. Otro inconveniente relevante del clustering de la liga sencilla es que podría ser difícil de interpretar en términos de particiones de los datos. Es una desventaja real del clustering de la liga sencilla porque estamos interesados en las particiones de los datos. Este método se tiende a desempeñar bien cuando hay grupos de forma elongada, conocidos como tipo cadena.
cluster_single <- hclust (bc_dist, method = 'single')
plot(cluster_single)
Cluster con el método single.
Clustering aglomerativo de la liga completa
El clustering de la liga completa también se conoce como clustering del vecino más lejano. Permite que que una muestra (o un grupo) se aglomere con otra muestra (o grupo) sólo a la distancia más lejana entre sí. Así, todos los miembros de ambos grupos están vinculados.
\[d_{AB}=\max_{i\in A, j\in B} (d_{ij}).\]
tabla de distancias
La agrupación de la liga completa evita el inconveniente de encadenar las muestras por el método de la liga simple. El encadenamiento completo tiende a encontrar muchos pequeños grupos separados. Este método se desempeña bien cuando los conglomerados son de forma circular.
cluster_complete <- hclust (bc_dist, method = 'complete')
plot(cluster_complete)
Cluster con el método complete.
Clustering aglomerativo de la liga promedio
La agrupación de liga promedio permite agrupar una muestra en un clúster en la media de las distancias entre esta muestra y todos los miembros del clúster. Los dos clusters se unen a la media de las distancias entre todos los miembros de un clúster y todos los miembros del otro.
\[d_{AB}=\frac{1}{n_An_B} \sum_{i\in A} \sum_{j\in B} (\ d_{ij}\ ).\]
tabla de distancias
El clúster resultante del dendrograma tiene un aspecto intermedio entre una agrupación de enlace simple y una completa.
cluster_average <- hclust (bc_dist, method = 'average')
plot(cluster_average)
Cluster con el método average.
Agrupación de varianza mínima de Ward
La agrupación de varianza mínima de Ward (también conocida como agrupación de Ward) presentada originalmente por Ward (1963) se basa en el criterio del modelo lineal de mínimos cuadrados. Este método método minimiza las sumas de las distancias al cuadrado (es decir, el error al cuadrado de ANOVA) entre las muestras. Básicamente, considera el análisis de conglomerados como un análisis de la varianza, en lugar de utilizar métricas de distancia o medidas de asociación. Este método es más apropiado para las variables cuantitativas, y no variables binarias. El clustering de Ward se implementa mediante la búsqueda del par de clusters en cada paso que conduce al mínimo incremento de la varianza total dentro del cluster después de la fusión. Este incremento es una distancia cuadrada ponderada entre los centros de los clusters. Aunque las distancias iniciales de los clusters se definen como la distancia euclidiana al cuadrado entre puntos en el clustering de Ward, otros métodos de distancia pueden producir resultados significativos.
tabla de distancias
cluster_ward <- hclust (bc_dist, method = 'ward.D2')
plot(cluster_ward)
Cluster con el método Ward.
# Correrlo en consola para visualizarlo mejor
par (mfrow = c(2,2))
plot(cluster_single)
plot(cluster_complete)
plot(cluster_average)
plot(cluster_ward)
Comparación de los clusters con los 4 métodos.
par (mfrow = c(1,1))La comparación entre estos cuatro dedrogramas muestra que la vinculación completa y Ward generan los mismos clusters de partición y tienen el mejor rendimiento en términos de partición de datos en la dirección de nuestro interés (las muestras 11, 12 13, 14 y 15 son de ratones Vdr\(-/-\), y 22, 23 y 24 de ratones de tipo salvaje).
Mismos gráficos usando el paquete factoextra
Heatmap de la matriz de distancias
set.seed(12345)
hc_single <- hclust(d=bc_dist, method = "single")
hc_average <- hclust(d=bc_dist, method = "average")
hc_complete <- hclust(d=bc_dist, method = "complete")
hc_ward <- hclust(d=bc_dist, method = "ward.D2")
con el método del promedio.
con el método Ward.
con el método de la liga más sencilla.
con el método de la liga completa.
Ordenación
El principal objetivo de la ordenación fue la “exploración” con la introducciónm de un análisis canonico de correspondencia (CCA), la ordenación ha ido más allá de meramente un análisis “exploratorio” y se ha convertido en una prueba de hipótesis.
La estructura de los datos de dos variables suele revelarse mediante un gráfico de dispersión de las muestras. Los datos multivariantes del microbioma son multidimensionales, generalmente tienen más de dos variables. El conjunto de datos del microbioma puede verse como una colección de muestras (sujetos) situadas en un espacio en el que cada variable o especie (o OTU/taxa) define una dimensión. Por lo tanto, hay tantas dimensiones como variables o especies (u OTU/taxa). Por ejemplo, en nuestro conjunto de datos fecales de ratón Vdr\(-/-\) hay 8 muestras y 248 géneros. Por tanto, la dimensión de los datos es 248. Para un conjunto de datos con n variables, el número de gráficos de dispersión que hay que dibujar sería \(n (n - 1)/2\). En nuestro caso, los 248 géneros necesitarán \((248*247)/2 = 30628\) gráficos de dispersión. Un número tan grande de gráficos de dispersión no es informativo para conocer la estructura de los datos; también es tedioso trabajar con ellos.
La ordenación trata principalmente de representar las relaciones entre las muestras y las especies (o OTUs/taxa) lo más fielmente posible en un espacio de baja dimensión (Gauch 1982a, b). Este objetivo es deseable, ya que los datos de la comunidad tienen múltiples dimensiones mezcladas con ruido, las dimensiones bajas pueden representar idealmente y de forma típica interpretaciones importantes e intuitivas de las relaciones especie (o OTUs/taxa)-ambiente. Para un conjunto de datos \(n\times p\) que contiene \(n\) sujetos y \(p\) variables, la estructura de los datos puede revisarse como \(n\) sujetos (representados como un grupo de puntos) en el espacio \(p-\)dimensional.
El objetivo principal de la ordenación es representar múltiples muestras (sujetos) en un número reducido de ejes ortogonales (es decir, independientes), donde el número total de ejes es menor o igual al número de muestras (sujetos). La importancia de los ejes de ordenación disminuye por orden. El primer eje de una ordenación explica la mayor parte de la variación en el conjunto de datos, seguido por el segundo eje, luego el tercero, y así sucesivamente. Los gráficos de ordenación son especialmente útiles para visualizar la similitud entre muestras (sujetos). Por ejemplo, en el contexto de la diversidad beta, las muestras que están más cerca en el espacio de ordenación tienen conjuntos de especies que son más similares entre sí que las muestras que están más alejadas en el espacio de ordenación.
Los métodos de ordenación incluyen ordenaciones restringidas y no restringidas, dependiendo de si los ejes de ordenación están limitados por factores ambientales (variables). Como su nombre indica, en la ordenación restringida, los ejes de ordenación están están restringidos por factores ambientales: las posiciones de las muestras en la ordenación están limitadas por las variables del entorno. En la ordenación no restringida, los ejes de ordenación ejes de ordenación no están limitados por factores ambientales. En otras palabras, la ordenación restringida utiliza directamente las variables ambientales en la construcción de la ordenación, mientras que en la ordenación no restringida los ejes no están limitados por factores ambientales; mientras que en la ordenación sin restricciones, las variables ambientales se introducen en los análisis post hoc.
Desde el punto de vista de la comprobación de hipótesis, el
análisis de ordenación sin restricciones per se es un método
principalmente descriptivo y no implica realmente la comprobación de
hipótesis en datos multivariados. Aunque implica la formulación de
hipótesis sobre la explicación de los ejes mediante variables
ambientales, como la función envfit() del paquete
vegan, la ordenación sin restricciones es sencilla. Analiza
una matriz de datos; su objetivo es revelar la estructura principal de
los datos en un gráfico a partir de un conjunto reducido de ejes
ortogonales. Por el contrario, la ordenación restringida es una
ordenación “impulsada por la hipótesis”: un método de comprobación de
hipótesis, que pone a prueba directamente la hipótesis sobre la
influencia de los factores ambientales en la composición de las especies
(o OTU/taxa).
La ordenación restringida está relacionada con los modelos lineales multivariantes, de esta manera con las variables “dependientes” (o la comunidad) en el lado izquierdo como respuestas, las “indenpendientes” (o restricciones) en el lado derecho como factores explicados. Así, la ordenación restringida no es simétrica.
En este apartado se tratan las ordenaciones no restringidas más comunes:
- Análisis de componentes principales (PCA).
- Análisis de coordenadas principales (PCoA).
- Escalamiento multidimensional no métrico (NMDS).
- Análisis de correspondencia (AC).
Y las ordenaciones restringidas:
- Análisis de redundancia redundancia (RDA).
- Análisis de correspondencia restringido (CCA).
- Análisis restringido de coordenadas principales (CAP)*.
Ordenaciones NO restringidas
Análisis de componentes principales (PCA)
En cuanto al paquete vegan, la variable
ambiental en nuestro caso es la condición genética (Vdr\(-/-\) y WT). Queremos saber si la
deficiencia del gen vdr interpreta la diversidad beta en la composición
del género. Realizamos un PCA para explorar si los cambios en
la composición de géneros de las comunidades (diversidades beta) son
causados por el factor genético. Con el PCA, las muestras se trazan en
base a las abundancias del género A en el eje 1 el género B en el eje 2,
el género C en el eje 3, y así sucesivamente hasta que se trazan \(n\) muestras en un espacio de muy alta
dimensión. Las \(n\) muestras crean
\((n - 1)\) número de componentes
principales (PC):
- PC1: es la primera línea recta que atraviesa el espacio creado por todas estas muestras.
- PC2: es la segunda línea, perpendicular a PC1.
y así la tercera PC3, hasta PC(n - 1). La importancia de las PC
disminuye por orden. PC1 es la PC más importante y
explica la mayor parte de las variaciones entre todas las muestras. La
PC2 explica la segunda mayor parte de las variaciones,
y PC3 explica la tercera más, y así sucesivamente la \((n - 1)\) PC explica la menor cantidad de
variaciones de las muestras. Varias funciones de R, como
prcomp() en el paquete preinstalado de estadísticas,
rda() en el paquete paquete vegan,
pca() en el paquete labdsv pueden utilizarse
para llevar a cabo el PCA.
En este caso, utilizamos la función rda(). Las dos
funciones de extensión: evplot() (Borcard et al. 2011) y
PCAsignificance() en el paquete BiodiversityR pueden
mejorar los gráficos. La función evplot() proporciona
métodos visuales para decidir la importancia de los ejes de ordenación
mediante el uso de el criterio de Keiser-Guttman y el modelo de
broken-stick. La función PCAsignificance() calcula el
modelo de broken-stick para los ejes PCA. Cuando se utiliza la función
rda() para realizar el PCA, al no especificar la matriz de
datos ambientales (es decir, la variable de grupo), la función realiza
un ordenación PCA. En el análisis de datos del microbioma, los recuentos
de abundancia absoluta no son apropiados porque los valores grandes
tendrán una influencia demasiado alta en el análisis. Por lo tanto,
necesitamos estandarizar los datos de lectura de abundancia
antes del análisis. En este caso, utilizamos la función
decostand() con el método total para estandarizar la
lectura.
# Llamamos a nuestra tabla
abund_table=read.csv(paste0(workingDir,"/data/VdrFecalGenusCounts.csv"),row.names=1,check.names=FALSE)
# Transpuesta de la tabla
abund_table<-t(abund_table)
#head(abund_table)
# Tabla de metadatos
meta_table <- data.frame(row.names=rownames(abund_table),
t(as.data.frame(strsplit(rownames(abund_table),"_"))))
# Agregamos el factor de grupo
meta_table$Group <- with(meta_table,ifelse(as.factor(X2)%in% c(11,12,13,14,15),
c("Vdr-/-"), c("WT")))
# Estandarizamos la abundancia de los datos
stand_abund_table <- decostand(abund_table, method = "total")
PCA <-rda(stand_abund_table)
PCA## Call: rda(X = stand_abund_table)
##
## Inertia Rank
## Total 0.04082
## Unconstrained 0.04082 7
## Inertia is variance
##
## Eigenvalues for unconstrained axes:
## PC1 PC2 PC3 PC4 PC5 PC6 PC7
## 0.020288 0.012792 0.003712 0.001730 0.001425 0.000774 0.000098Con el siguiente código podemos pedir los valores de los eigenvalores.
ev <- PCA$CA$eigEn lenguaje vegan, “Inertia” es el
término general de “variación” en los datos. El total
de variación de todo el conjunto de datos es de 0.0408196 en este caso,
y el primer eje explica el 49.7017233% de la variación total
(0,02029/0,0408 = 0,4973). La variación total es la
suma de las variaciones de cada género en la matriz analizada. Los
siguientes códigos R comprueban la variación total:
# Variación total
sum(apply(stand_abund_table, 2, var))## [1] 0.04081957Otra forma de obtener este análisis es llamando la función
summary.
head(summary(PCA))##
## Call:
## rda(X = stand_abund_table)
##
## Partitioning of variance:
## Inertia Proportion
## Total 0.04082 1
## Unconstrained 0.04082 1
##
## Eigenvalues, and their contribution to the variance
##
## Importance of components:
## PC1 PC2 PC3 PC4 PC5 PC6
## Eigenvalue 0.02029 0.01279 0.003712 0.00173 0.001425 0.0007736
## Proportion Explained 0.49702 0.31338 0.090938 0.04239 0.034921 0.0189506
## Cumulative Proportion 0.49702 0.81039 0.901331 0.94372 0.978643 0.9975932
## PC7
## Eigenvalue 9.824e-05
## Proportion Explained 2.407e-03
## Cumulative Proportion 1.000e+00
##
## Scaling 2 for species and site scores
## * Species are scaled proportional to eigenvalues
## * Sites are unscaled: weighted dispersion equal on all dimensions
## * General scaling constant of scores: 0.731125
##
##
## Species scores
##
## PC1 PC2 PC3 PC4
## Tannerella -0.1460060 0.059004 -0.1581134 0.0263299
## Lactococcus 0.3347414 0.268815 0.0199084 -0.0094492
## Lactobacillus 0.2884642 -0.241439 -0.0576255 0.0148015
## Lactobacillus::Lactococcus 0.0047630 0.002077 -0.0002720 -0.0008367
## Parasutterella 0.0002188 -0.001427 -0.0000296 -0.0008465
## Helicobacter -0.0390189 0.005225 -0.0389887 -0.0101868
## ....
## PC5 PC6
## Tannerella -0.0074351 -0.0108844
## Lactococcus 0.0133844 -0.0002257
## Lactobacillus 0.0420286 0.0006774
## Lactobacillus::Lactococcus 0.0037431 -0.0004606
## Parasutterella -0.0009601 0.0001548
## Helicobacter -0.0106939 -0.0092049
## ....
##
##
## Site scores (weighted sums of species scores)
##
## PC1 PC2 PC3 PC4 PC5 PC6
## 5_15_drySt-28F -0.35858 0.04972 -0.443749 -0.2395 -0.12631 -0.21700
## 1_11_drySt-28F -0.01508 0.25950 0.308436 -0.4460 0.28194 0.14170
## 2_12_drySt-28F -0.28343 -0.36080 0.318538 0.1051 0.10647 -0.28342
## 3_13_drySt-28F 0.02043 0.37303 -0.007325 0.1292 -0.26490 0.07505
## 4_14_drySt-28F -0.24988 -0.01016 0.007908 0.3703 0.03704 0.44259
## 7_22_drySt-28F 0.31464 -0.02094 -0.313787 0.1650 0.47544 -0.09180
## 8_23_drySt-28F 0.28787 -0.42150 -0.065406 -0.2454 -0.26694 0.26212
## 9_24_drySt-28F 0.28403 0.13115 0.195385 0.1613 -0.24274 -0.32924
Explicación PCA
También podemos obtener cuales variables (species)
tienen la mayor correlación absoluta al primer y segundo eje:
head(sort(abs(PCA$CA$v[,1]),decreasing = TRUE))## Lactococcus Lactobacillus Tannerella Clostridium Bacteroides
## 0.6494296 0.5596474 0.2832653 0.2751724 0.2046105
## Allobaculum
## 0.1487113head(sort(abs(PCA$CA$v[,2]),decreasing = TRUE))## Lactococcus Lactobacillus Allobaculum Akkermansia Bacteroides
## 0.6567934 0.5899062 0.3501612 0.1680896 0.1457728
## Tannerella
## 0.1441641A continuación, vamos a dibujar los diagramas utilizando la función
biplot(). La opción de visualización “species” es la
etiqueta del paquete vegan para OTUs/taxa. Por defecto es
“sites” (etiqueta para muestras).
biplot(PCA, display = 'species')
Biplot de las primeras dos componentes principales de los datos de materia fecal de ratones Vdr\(-/-\)
Otros argumentos que uno puede incluir en el biplot son por ejemplo
mostrar también los sites y agregar etiquetas a los datos,
(Ver vegan
pag. 37)
#Mismo biplot con etiquetas
biplot(PCA, display = c('sites','species'), type=c("text","points"))
Biplot con etiquetas
Los diagramas anteriores trazados por biplot() sólo
dibujan flechas para el género, lo que no es informativo. El diagrama
más informativo es utilizar la función ordiplot() para
dibujar tanto el género como las puntuaciones de la muestra como
centroides, como se muestra a continuación:
ordiplot(PCA, display = "sites", type = "text")
Ordiplot de las primeras dos componentes principales de los datos de materia fecal de ratones Vdr\(-/-\) con etiquetas
En los argumentos anteriores, type="text" o
t añade etiquetas de texto a la figura (la configuración
por defecto sólo añade puntos). Como alternativa, podemos utilizar la
función cleanplot.pca() escrita por François Gillet &
Daniel Borcard para intentar dibujar los resultados del PCA en dos
diagramas y diferenciarlos por su escala. Esta función sirve para
dibujar dos biplots (escalado 1 y escalado 2) a partir de un objeto de
clase “rda” en PCA o RDA resultado de la función rda() de
vegan. Primero ejecutamos la función cleanplot.pca() como
se indica a continuación.
"cleanplot.pca" <- function(res.pca, ax1=1, ax2=2, point=FALSE,
ahead=0.07, cex=0.5)
{
# A function to draw two biplots (scaling 1 and scaling 2) from an object
# of class "rda" (PCA or RDA result from vegan's rda() function)
#
# License: GPL-2
# Authors: Francois Gillet & Daniel Borcard, 24 August 2012
require("vegan")
par(mfrow=c(1,2))
p <- length(res.pca$CA$eig)
# Scaling 1: "species" scores scaled to relative eigenvalues
sit.sc1 <- scores(res.pca, display="wa", scaling=1, choices=c(1:p))
spe.sc1 <- scores(res.pca, display="sp", scaling=1, choices=c(1:p))
plot(res.pca, choices=c(ax1, ax2), display=c("wa", "sp"), type="n",
main="PCA - scaling 1", scaling=1)
if (point)
{
points(sit.sc1[,ax1], sit.sc1[,ax2], pch=20)
text(res.pca, display="wa", choices=c(ax1, ax2), cex=cex, pos=1, scaling=1)
}
else
{
text(res.pca, display="wa", choices=c(ax1, ax2), cex=cex, scaling=1)
}
text(res.pca, display="sp", choices=c(ax1, ax2), cex=cex, pos=1, col="blue", scaling=1)
arrows(0, 0, spe.sc1[,ax1], spe.sc1[,ax2], length=ahead, angle=20, col="red")
pcacircle(res.pca)
# Scaling 2: site scores scaled to relative eigenvalues
sit.sc2 <- scores(res.pca, display="wa", choices=c(1:p))
spe.sc2 <- scores(res.pca, display="sp", choices=c(1:p))
plot(res.pca, choices=c(ax1,ax2), display=c("wa","sp"), type="n",
main="PCA - scaling 2")
if (point) {
points(sit.sc2[,ax1], sit.sc2[,ax2], pch=20)
text(res.pca, display="wa", choices=c(ax1 ,ax2), cex=cex, pos=1)
}
else
{
text(res.pca, display="wa", choices=c(ax1, ax2), cex=cex)
}
text(res.pca, display="sp", choices=c(ax1, ax2), cex=cex, pos=1, col="blue")
arrows(0, 0, spe.sc2[,ax1], spe.sc2[,ax2], length=ahead, angle=20, col="red")
}
"pcacircle" <- function (pca)
{
# Draws a circle of equilibrium contribution on a PCA plot
# generated from a vegan analysis.
# vegan uses special constants for its outputs, hence
# the 'const' value below.
eigenv <- pca$CA$eig
p <- length(eigenv)
n <- nrow(pca$CA$u)
tot <- sum(eigenv)
const <- ((n - 1) * tot)^0.25
radius <- (2/p)^0.5
radius <- radius * const
symbols(0, 0, circles=radius, inches=FALSE, add=TRUE, fg=2)
}cleanplot.pca(PCA)
PCA plots de las primeras dos componentes principales de los datos de materia fecal de ratones Vdr\(-/-\) usando la función cleant.pca()
El “escalado” es la forma en que los resultados de la
ordenación se proyectan en el espacio reducido para su visualización
gráfica (Borcard et al. 2011). La función cleanplot.pca()
genera dos escalas de PCA. El escalado PCA 1 en la figura de la
izquierda es el biplot de distancia que se centra en las distancias
entre las muestras. Si está interesado principalmente en interpretar las
relaciones entre las muestras, elija el escalado 1. Las características
esenciales de Scaling 1 son: Los vectores propios se escalan a la unidad
de longitud y las distancias entre las muestras (sujetos) en el biplot
son aproximaciones de sus distancias euclidianas en el espacio
multidimensional. Sin embargo, los ángulos entre los vectores de las
variables (bacterias en este caso) no tienen sentido. El círculo se
llama círculo de contribución de equilibrio, que representa la
contribución de equilibrio de las variables. Para una combinación dada
de ejes, las variables con vectores más largos que el radio del
círculo pueden ser interpretadas con confianza como las bacterias
más importantes, mientras que las variables con vectores más
cortos que el radio del círculo de contribución de equilibrio
contribuyen poco a un espacio reducido dado.
El PCA 2 de la figura de la derecha es un biplot de correlación. Traza la correlación entre variables (bacterias en este caso). Si su interés principal se centra en las relaciones entre variables, elija el escalamiento 2. Las características esenciales del escalamiento 2 son: Cada vector propio se escala a la raíz cuadrada de su valor propio; la longitud del vector se aproxima a la desviación estándar de las variables (bacterias); los ángulos entre variables (bacterias en este caso) reflejan sus correlaciones: el coseno del ángulo aproxima la correlación entre las variables (bacterias); sin embargo, las distancias entre muestras (sujetos) en el biplot no son aproximaciones de sus distancias euclidianas en el espacio multidimensional.
Cálculo del PC1
El análisis de componentes principales tiene como objetivo reducir la dimensión y conservar en lo posible la estructura de estos, la cual no depende de los ejes utilizados para las coordenadas. Dadas \(n\) observaciones de \(p\) variables, se analiza si es posible representar adecuadamente la información con un número menor de variables construidas como combinaciones lineales de las originales. El primer componente principal es la combinación lineal de las variables originales que tienen varianza máxima.
Supongamos que tenemos datos en un espacio de \(p\) variables en \(n\) elementos de una población en una matriz \(X\) de tamaño \(n\times p\), donde las columnas son las variables y las filas los elementos. Vamos a suponer que hemos restado a cada variable su media, de tal forma que las variables de la matriz \(X\) tienen media cero y su matriz de covarianzas está dada por \(\frac{1}{n}X'X\). Los valores del primer componente de los \(n\) individuos se representará por un vector \(\mathbb{z}_1\) que está dado por: \[\mathbb{z}_1=Xa_{1}.\]
Donde \(a_1=(a_{11},\ a_{12},...,\ a_{1p})\). Como las variables originales se supuso tenian media cero, entonces también \(\mathbb{z}_1\) tendrá media nula y su varianza sera: \[\frac{1}{n}\ \mathbb{z}_{1}'\ \mathbb{z}_1 = \frac{1}{n}\ a_{1}X'Xa_1\ =\ a_1' S a_1.\] donde \(S\) es la matriz de varianzas y covarianzas de las observaciones. Si se elige un \(a_1\) que tenga una norma muy grande, entonces la varianza de \(\mathbb{z}_1\) podría resultad muy grande también. Entonces es necesario imponer condiciones sobre \(a_1\) para acotar el tamaño de la varianza de \(\mathbb{z}_1\) y no hacerla muy grande de manera arbitraria. La condición que se impone es: \[||a_1||=a_1'\cdot a_1=1.\]
Esta restricción se introduce mediante el multiplicador de Lagrange:
\[M\ =\ a_1' S a_1\ -\ \lambda(a_1'a_1\ -\ 1). \]
Después, se maximizará esta expresión derivando con respecto a los componentes de \(a_1\) e igualando a \(0\):
\[ \frac{\partial M}{\partial a_1} =\ 2Sa_1 -2\lambda a_1 =0 \]
cuya solución es:
\[ Sa_1=\lambda a_1 \]
entonces \(a_1\) es un vector propio de la matriz \(S\) y \(\lambda\) será su correspondiente valor propio.
Para determinar que valor propio de \(S\) es la solución de la ecuación anterior, multiplicamos por \(a_1\) a la izquierda:
\[a_1'Sa_1=\lambda a_1'a_1 =\lambda \]
entonces, \(\lambda\) es la varianza de \(\mathbb{z}_1\). Como esta es la cantidad que queríamos maximizar, entonces \(\lambda\) será el valor propio mayor de la matriz \(S\) y su vector asociado \(a_1\) serán los coeficientes de cada variable de la primera componente.
Análisis de coordenadas principales (PCoA)
El PCoA también se denomina escalamiento multidimensional
métrico. PCoA es una técnica de ordenación que permite al usuario
elegir prácticamente cualquier métrica de distancia (por ejemplo,
Jaccard, Bray-Curtis, Euclides, etc.). Al
igual que PCA, PCoA utiliza los valores propios para medir la
importancia de un conjunto de ejes ortogonales devueltos. La
dimensionalidad de la matriz se reduce determinando cada vector propio y
cada valor propio. Las coordenadas principales se obtienen escalando
cada vector propio. El PCoA cuando se calcula sobre distancias
euclidianas entre las muestras produce los mismos resultados que el PCA
calculado sobre matriz de covarianza del mismo conjunto de datos (si se
utiliza el escalado 1). Las funciones de R, incluyendo
cmdscale() en el paquete vegan y
pcoa() en el paquete ape pueden realizar el
PCoA. Con vegan, los datos de entrada pueden ser calculados por la
función vegdist() (por defecto es la disimilitud de
Bray-Curtis), y el diagrama de ordenación puede ser dibujado
con la función ordiplot(). El diagrama de ordenación
también puede ser con la función biplot.pcoa() del paquete
ape. En este caso, utilizamos la función
cmdscale() y los mismos datos fecales de ratones Vdr\(-/-\) para realizar un PCoA. Esta función
necesita una matriz de semejanza como datos de entrada. Primero, vamos a
calcular la disimilitud de Bray-Curtis utilizando la función
vegdist().
options(digits=4)
options(width=78, digits=4)
library("vegan")
bc_dist <-vegdist(abund_table, "bray")| 5_15_ | 1_11_ | 2_12_ | 3_13_ | 4_14_ | 7_22_ | 8_23_ | 9_24_dry | |
|---|---|---|---|---|---|---|---|---|
| 5_15_drySt-28F | 0.0000 | 0.5873 | 0.5293 | 0.5786 | 0.2765 | 0.4756 | 0.6487 | 0.5984 |
| 1_11_drySt-28F | 0.5873 | 0.0000 | 0.3742 | 0.2041 | 0.4896 | 0.5208 | 0.4084 | 0.2921 |
| 2_12_drySt-28F | 0.5293 | 0.3742 | 0.0000 | 0.4861 | 0.4114 | 0.4955 | 0.5119 | 0.4368 |
| 3_13_drySt-28F | 0.5786 | 0.2041 | 0.4861 | 0.0000 | 0.4645 | 0.5436 | 0.4217 | 0.2362 |
| 4_14_drySt-28F | 0.2765 | 0.4896 | 0.4114 | 0.4645 | 0.0000 | 0.4172 | 0.6220 | 0.5210 |
| 7_22_drySt-28F | 0.4756 | 0.5208 | 0.4955 | 0.5436 | 0.4172 | 0.0000 | 0.5356 | 0.4359 |
| 8_23_drySt-28F | 0.6487 | 0.4084 | 0.5119 | 0.4217 | 0.6220 | 0.5356 | 0.0000 | 0.3105 |
| 9_24_drySt-28F | 0.5984 | 0.2921 | 0.4368 | 0.2362 | 0.5210 | 0.4359 | 0.3105 | 0.0000 |
A continuación, vamos a establecer explícitamente k = 2
(los valores por defecto para el número de dimensiones a devolver) y
eig = TRUE (que guarda los valores propios).
PCoA <- cmdscale(bc_dist, eig = TRUE, k = 2)
PCoA## $points
## [,1] [,2]
## 5_15_drySt-28F 0.35060 0.007035
## 1_11_drySt-28F -0.17119 0.137244
## 2_12_drySt-28F 0.02928 0.115931
## 3_13_drySt-28F -0.17013 0.126964
## 4_14_drySt-28F 0.27190 0.088394
## 7_22_drySt-28F 0.13943 -0.258936
## 8_23_drySt-28F -0.25091 -0.157891
## 9_24_drySt-28F -0.19898 -0.058740
##
## $eig
## [1] 3.779e-01 1.517e-01 1.170e-01 8.855e-02 2.771e-02 2.167e-02
## [7] -2.515e-17 -4.493e-03
##
## $x
## NULL
##
## $ac
## [1] 0
##
## $GOF
## [1] 0.6713 0.6751La función cmdscale() produce una lista de salida. Los
primeros puntos de salida contienen las coordenadas de
cada muestra en cada dimensión reducida. La segunda salida
eig contiene los valores propios. Las últimas tres
salidas pertenecen a otras opciones del análisis que no cubriremos aquí.
El siguiente fragmento de códigos R se utiliza para examinar el
porcentaje de variación en el conjunto de datos que se explica por los
dos primeros ejes del PCoA.
explainedvar1 <- round(PCoA$eig[1] / sum(PCoA$eig), 2) * 100
explainedvar1## [1] 48explainedvar2 <- round(PCoA$eig[2] / sum(PCoA$eig), 2) * 100
explainedvar2## [1] 19sum_eig <- sum(explainedvar1, explainedvar2)
sum_eig## [1] 67El primer eje explica el 48% de las variaciones de los datos, y el segundo el 19%. Por lo tanto, una gran cantidad de variaciones de los datos (un total del 67%) ha sido explicada por estos dos ejes.
Existen dos criterios para evaluar si los primeros ejes del PCoA captan o no una cantidad desproporcionadamente grande de la variación total explicada: (1) el criterio de Kaiser-Guttman y (2) el modelo broken-stick. El criterio de Kaiser-Guttman establece que “los valores propios asociados a los primeros ejes deben ser mayores que la media de todos los valores propios”. Según el criterio del modelo broken-stick, los valores propios asociados a los primeros ejes se comparan con las expectativas del modelo broken-stick. El modelo de broken-stick supone que la suma total de los valores propios disminuye secuencialmente con los ejes PCoA ordenados. Evaluamos el rendimiento de PCoA utilizando estos dos criterios con los siguientes gráficos:
# Correr todo junto
# Definir parámetros
par(mar = c(5, 5, 1, 2) + 0.1)
# Graficar eigenvalores
plot(PCoA$eig, xlab = "PCoA", ylab = "Eigenvalue",
las = 1, cex.lab = 1.5, pch = 16)
#Añadir expextación del criterio de Kaiser-Guttman y el modelo Broken Stick
abline(h = mean(PCoA$eig), lty = 2, lwd = 2, col = "blue")
b_stick <- bstick(8, sum(PCoA$eig))
lines(1:8, b_stick, type = "l", lty = 4, lwd = 2, col = "red")
# Añadir legendas
legend("topright", legend = c("Avg Eigenvalue", "Broken-Stick"),
lty = c(2, 4), bty = "n", col = c("blue", "red"))
PCoA con el criterio de Kaiser-Guttman y el modelo de broken-stick
La figura anterior muestra que los valores propios asociados a los tres primeros ejes son mayores que la media de todos los valores propios y son mayores que las expectativas del el modelo Broken-Stick. Después de evaluar el resultado del PCoA, crearemos un de ordenación para los dos ejes del PCoA.
# Correr todo junto
# Definir parámetros
par(mar = c(5, 5, 1, 2) + 0.1)
# Inicializar el plot
plot(PCoA$points[ ,1], PCoA$points[ ,2], ylim = c(-0.5, 0.5),
xlab = paste("PCoA 1 (", explainedvar1, "%)", sep = ""),
ylab = paste("PCoA 2 (", explainedvar2, "%)", sep = ""),
pch = 5, cex = 1.0, type = "n", cex.lab = 1.0, cex.axis = 1.2, axes = FALSE)
# Añadir ejes
axis(side = 1, labels = T, lwd.ticks = 2, cex.axis = 1.2, las = 1)
axis(side = 2, labels = T, lwd.ticks = 2, cex.axis = 1.2, las = 1)
abline(h = 0, v = 0, lty = 3)
box(lwd = 2)
# Añadir puntos y etiquetas
points(PCoA$points[ ,1], PCoA$points[ ,2],
pch = 19, cex = 3, bg = "blue", col = "blue")
text(PCoA$points[ ,1], PCoA$points[ ,2],
labels = row.names(PCoA$points))
Ordination plot para los dos ejes de PCoA
Los gráficos de ordenación básicos nos permiten ver cómo se separan
las muestras entre sí. En nuestro ejemplo, las muestras se separan a lo
largo de los ejes PCoA debido a la variación en la abundancia de los
diferentes géneros de ratones. Una pregunta lógica de seguimiento es
preguntar qué géneros del conjunto de datos están impulsando la
divergencia observada entre los puntos. ¿Podemos identificar y
visualizar estos géneros influyentes en el PCoA? Podemos obtener esta
información utilizando la función add.spec.scores() del
paquete BiodiversityR. En primer lugar, la abundancia relativa se
calcula como sigue:
# Correr todo junto
# Definir parámetros del plot
par(mar = c(5, 5, 1, 2) + 0.1)
# Inicializar el plot
plot(PCoA$points[ ,1], PCoA$points[ ,2], ylim = c(-0.5, 0.5),
xlab = paste("PCoA 1 (", explainedvar1, "%)", sep = ""),
ylab = paste("PCoA 2 (", explainedvar2, "%)", sep = ""),
pch = 5, cex = 1.0, type = "n", cex.lab = 1.0, cex.axis = 1.2, axes = FALSE)
# Añadir ejes
axis(side = 1, labels = T, lwd.ticks = 2, cex.axis = 1.2, las = 1)
axis(side = 2, labels = T, lwd.ticks = 2, cex.axis = 1.2, las = 1)
abline(h = 0, v = 0, lty = 3)
box(lwd = 2)
# Añadir puntos y etiquetas
points(PCoA$points[ ,1], PCoA$points[ ,2],
pch = 19, cex = 3, bg = "blue", col = "blue")
text(PCoA$points[ ,1], PCoA$points[ ,2],
labels = row.names(PCoA$points))
fecalREL <- abund_table
for(i in 1:nrow(abund_table)){
fecalREL[i, ] = abund_table[i, ] / sum(abund_table[i, ])
}
#install.packages("BiodiversityR")
library("pbkrtest")
library("BiodiversityR")
# Calcular y añadir el score de las especies
PCoA <- add.spec.scores(PCoA,fecalREL,method = "pcoa.scores",Rscale=TRUE,scaling=1, multi=1)
text(PCoA$cproj[ ,1], PCoA $cproj[ ,2],
labels = row.names(PCoA$cproj),cex=0.5, col = "blue")
Ordination plot para los dos ejes de PCoA con influencia de genero
La función add.spec.scores() también puede utilizarse
para determinar la correlación de cada taxón a lo largo de los ejes
PCoA. Este es un enfoque más efectivamente cuantitativo de identificar
los taxones influyentes. Además, para identificar y obtener el taxón más
importante, es necesario definir un coeficiente de correlación de corte,
en este caso, vamos a tomar \(0.7\).
genus_corr <- add.spec.scores(PCoA, fecalREL, method = "cor.scores")$cproj
corrcut <- 0.7 # definir el cutoff
import_genus <- genus_corr[abs(genus_corr[, 1]) >= corrcut | abs(genus_corr[, 2]) >= corrcut, ]Los 12 géneros importantes con correlación mayor o igual a 0,7 a lo largo de los ejes PCoA se imprimen a continuación:
import_genus[complete.cases(import_genus),] %>%
kbl(booktabs = TRUE, align = "c",
caption = "Géneros más importantes con correlación mayor o igual a 0.7" )%>%
kable_styling(position = "center",
latex_options = c("hold_position", "striped"),
font_size = 12)| Dim1 | Dim2 | |
|---|---|---|
| Tannerella | 0.8205 | 0.1685 |
| Lactobacillus | -0.3064 | -0.8665 |
| Helicobacter | 0.7488 | 0.1176 |
| Paraprevotella | 0.7463 | -0.0858 |
| Bacillus | 0.1980 | -0.7727 |
| Pedobacter | -0.2461 | 0.8703 |
| Limibacter | -0.0824 | 0.7036 |
| Mycoplasma | 0.7227 | 0.2295 |
| Slackia | 0.0977 | -0.7577 |
| Fluviicola | 0.2425 | -0.7107 |
| Caldicellulosiruptor | 0.2425 | -0.7107 |
| Anaerophaga | 0.2425 | -0.7107 |
Por último, utilizamos la función envfit() del paquete
vegan para realizar una prueba de permutación para las
abundancias generales a través de los ejes en estas correlaciones.
# Test de permutación para abundancia de especies en los ejes
envfit(PCoA, fecalREL, perm = 999)##
## ***VECTORS
##
## Dim1 Dim2
## Tannerella 0.951 0.308
## Lactococcus -0.769 -0.639
## Lactobacillus -0.219 -0.976
## Lactobacillus::Lactococcus -0.428 -0.904
## Parasutterella -0.754 -0.657
## Helicobacter 0.971 0.241
## Prevotella 0.842 0.540
## Bacteroides 0.262 0.965
## Barnesiella 0.979 0.204
## Odoribacter -0.537 -0.843
## Eubacterium 0.238 0.971
## Allobaculum -0.676 -0.737
## Roseburia -0.044 0.999
## Clostridium 0.280 0.960
## Porphyromonas 0.359 -0.933
## Butyrivibrio -0.006 1.000
## Ruminococcus 0.557 0.830
## Acholeplasma 0.327 0.945
## Alistipes 0.071 0.997
## Clostridium::Coprococcus 0.000 0.000
## Eubacterium (Erysipelotrichaceae)::Eubacterium -0.820 0.573
## Hydrogenoanaerobacterium 0.037 -0.999
## Paraprevotella 0.984 -0.179
## Blautia -0.087 -0.996
## Adlercreutzia::Asaccharobacter -0.140 0.990
## Coprococcus -0.011 1.000
## Parabacteroides 0.978 0.211
## Eubacterium (Erysipelotrichaceae) -0.770 -0.639
## Oscillibacter 0.273 0.962
## Acinetobacter 0.999 0.050
## Alkalitalea -0.622 -0.783
## Sporobacter 0.000 0.000
## Oscillospira 0.481 0.876
## Blautia::Clostridium 0.107 -0.994
## Planctomyces -0.474 0.881
## Akkermansia -0.587 0.809
## Ruminofilibacter -0.538 -0.843
## Pseudoflavonifractor 0.639 0.769
## Butyricimonas 0.293 -0.956
## Desulfovibrio 0.971 0.237
## Anaerotruncus -0.448 0.894
## Alistipes::Bacteroides 0.009 1.000
## Turicibacter 0.293 0.956
## Mucispirillum 0.999 0.050
## Lachnospira 0.494 0.869
## Catabacter 0.000 0.000
## Desulfosporosinus 0.000 0.000
## Bacillus 0.160 -0.987
## Sporanaerobacter -0.585 -0.811
## Streptomyces -0.677 0.736
## Butyrivibrio::Clostridium 0.922 -0.388
## Candidatus Arthromitus -0.185 -0.983
## Enterorhabdus 0.366 0.931
## Kurthia 0.000 0.000
## Eggerthella 0.974 0.227
## Actinocorallia 0.000 0.000
## Caldanaerocella 0.999 0.050
## Dorea 0.391 -0.921
## Streptococcus -0.370 0.929
## Adlercreutzia -0.791 0.612
## Paraeggerthella 0.962 0.273
## Lachnobacterium -0.538 -0.843
## TM7 (genus) 0.101 0.995
## Clostridium::Anaerostipes 0.999 0.050
## Erysipelothrix -0.208 -0.978
## Formosa 0.000 0.000
## Rikenella -0.521 -0.854
## Dysgonomonas -0.456 0.890
## Desulfitobacterium 0.999 0.050
## Oribacterium -0.474 0.881
## Syntrophococcus 0.000 0.000
## Ruminococcus::Clostridium -0.538 -0.843
## Rhodospirillum 0.999 0.050
## Pedobacter -0.176 0.984
## Acetivibrio 0.996 -0.085
## Clostridium::Eubacterium -0.997 0.078
## Limibacter -0.074 0.997
## Clostridium::Ruminococcus 0.000 0.000
## Coprobacillus 0.990 0.142
## Cytophaga -0.403 -0.915
## Denitrobacterium -0.382 -0.924
## Mycoplasma 0.894 0.448
## Roseburia::Clostridium -0.448 0.894
## Sporobacterium 0.583 0.813
## Eubacterium::Ruminococcus 0.000 0.000
## Pseudobutyrivibrio 0.101 0.995
## Robinsoniella 0.000 0.000
## Brevibacterium -0.192 0.981
## Blautia::Ruminococcus 0.961 -0.278
## Pseudomonas -0.405 0.914
## Clostridium::Roseburia 0.000 0.000
## Desulfotomaculum -0.991 0.134
## Clostridium (Erysipelotrichaceae) 0.764 0.645
## Olsenella -0.714 0.700
## Azospirillum 0.000 0.000
## Oxobacter 0.000 0.000
## Asaccharobacter::Adlercreutzia 0.000 0.000
## Desulfocurvus 0.101 0.995
## Anaerostipes 0.101 0.995
## Clostridium::Ruminococcus::Desulfotomaculum::Escherichia 0.000 0.000
## Halanaerobacter -0.999 0.034
## Slackia 0.081 -0.997
## Anaplasma 0.000 0.000
## Syntrophomonas 0.000 0.000
## Clostridium::Blautia 0.000 0.000
## Anaerosporobacter::Clostridium 0.000 0.000
## Clostridium::Bacteroides 0.000 0.000
## Blautia::Lactonifactor 0.000 0.000
## Parabacteroides::Bacteroides -0.474 0.881
## Enterorhabdus::Adlercreutzia 0.000 0.000
## Flavobacterium -0.266 0.964
## Clostridium::Desulfotomaculum 0.000 0.000
## Parasporobacterium 0.999 0.050
## Coprococcus::Clostridium 0.000 0.000
## Fusobacterium 0.854 0.521
## Ruminococcus::Escherichia -0.448 0.894
## Clostridium::Ruminococcus::Coprococcus 0.000 0.000
## Lactonifactor 0.000 0.000
## Haemophilus 0.309 0.951
## Sporichthya -0.849 0.529
## Cytophaga::Flavobacterium -0.768 -0.640
## Clostridium::Dorea 0.000 0.000
## Ruminococcus::Blautia 0.000 0.000
## Frigoribacterium -0.474 0.881
## Paenibacillus -0.738 -0.675
## Johnsonella 0.101 0.995
## Pseudoflavonifractor::Clostridium 0.000 0.000
## Adlercreutzia::Enterorhabdus::Asaccharobacter 0.000 0.000
## Nocardiopsis -0.474 0.881
## Pedobacter::Pseudomonas 0.000 0.000
## Flexibacter -0.688 -0.726
## Catabacter::Ruminococcus 0.000 0.000
## Butyricimonas::Bacteroides -0.474 0.881
## Staphylococcus -0.253 0.968
## Alkalitalea::Prevotella 0.000 0.000
## Lachnospira::Anaerostipes 0.000 0.000
## Flavobacterium::Cytophaga -0.538 -0.843
## Gelidibacter 0.817 0.577
## Treponema -0.474 0.881
## Pontibacter 0.000 0.000
## Desulfuromonas 0.000 0.000
## Butyricicoccus 0.000 0.000
## Clostridium::Butyrivibrio -0.448 0.894
## Coprobacillus::Clostridium 0.000 0.000
## Porphyromonas::Prevotella 0.000 0.000
## Effluviibacter -0.580 0.815
## Caminicella 0.101 0.995
## Prevotella::Bacteroides 0.000 0.000
## Pseudoalteromonas 0.000 0.000
## Butyrivibrio::Blautia 0.000 0.000
## Lachnospira::Clostridium -0.806 -0.593
## Alicyclobacillus 0.000 0.000
## Lachnospira::Robinsoniella 0.000 0.000
## Subdoligranulum 0.000 0.000
## Anaerofilum 0.000 0.000
## Sporosarcina 0.000 0.000
## Alkalitalea::Sphingobacterium 0.000 0.000
## Kopriimonas 0.777 0.629
## Olivibacter 0.000 0.000
## Lachnobacterium::Coprococcus 0.000 0.000
## Fluviicola 0.211 -0.977
## Peptococcus 0.000 0.000
## Ruminococcus::Roseburia 0.000 0.000
## Marinilabilia 0.000 0.000
## Catonella 0.000 0.000
## Sphingobium 0.101 0.995
## Olsenella::Streptomyces 0.000 0.000
## Terasakiella 0.000 0.000
## Roseospirillum -0.152 -0.988
## Clostridium::Ruminococcus::Blautia 0.000 0.000
## Lachnospira::Roseburia 0.000 0.000
## Thalassospira 0.000 0.000
## Eubacterium::Clostridium 0.101 0.995
## Allobaculum::Eubacterium 0.000 0.000
## Blautia::Lachnospira 0.000 0.000
## Cryptobacterium 0.000 0.000
## Atopobium 0.000 0.000
## Eubacterium::Lachnospira 0.000 0.000
## Bacteroides::Alistipes 0.000 0.000
## Faecalibacterium -0.192 0.981
## Lachnospira::Coprococcus 0.000 0.000
## Microbacterium 0.000 0.000
## Zhangella::Stella 0.777 0.629
## Paludibacter 0.000 0.000
## Butyrivibrio::Ruminococcus 0.000 0.000
## Bacteroides::Lactobacillus 0.000 0.000
## Prevotella::Flavobacterium -0.806 -0.593
## Ruminococcus::Dorea 0.000 0.000
## Slackia::Asaccharobacter::Adlercreutzia 0.000 0.000
## Caldicellulosiruptor 0.211 -0.977
## Kordia 0.000 0.000
## Bacteroides::Lactobacillus::Lachnospira -0.448 0.894
## Sphingobacterium 0.000 0.000
## Anaeroplasma -0.448 0.894
## Atopostipes 0.000 0.000
## Enterococcus 0.000 0.000
## Insolitispirillum 0.000 0.000
## Clostridium::Acetivibrio 0.000 0.000
## Plantactinospora 0.000 0.000
## Stappia 0.000 0.000
## Anaplasma::Clostridium 0.000 0.000
## Clostridium (Erysipelotrichaceae)::Clostridium 0.000 0.000
## Algoriphagus 0.000 0.000
## Dorea::Ruminococcus 0.000 0.000
## Roseburia::Lachnospira 0.000 0.000
## Rhizobium 0.000 0.000
## Anaerofustis 0.000 0.000
## Echinicola 0.000 0.000
## Anaerostipes::Clostridium 0.000 0.000
## Lachnospira::Blautia 0.000 0.000
## Aestuariimicrobium 0.000 0.000
## Gelidibacter::Flavobacterium 0.000 0.000
## Helicobacter::Flexispira 0.101 0.995
## Eubacterium (Erysipelotrichaceae)::Paenibacillus::Eubacterium 0.000 0.000
## Pelospora 0.000 0.000
## Stella 0.000 0.000
## Methylobacterium 0.000 0.000
## Rickettsia -0.448 0.894
## Porphyromonas::Flavobacterium 0.000 0.000
## Adlercreutzia::Enterorhabdus 0.000 0.000
## Ruminococcus::Escherichia::Parabacteroides 0.000 0.000
## Limibacter::Bacteroides 0.000 0.000
## Paraprevotella::Prevotella 0.000 0.000
## Janthinobacterium::Zoogloea::Duganella 0.000 0.000
## Flavobacterium::Gelidibacter 0.000 0.000
## Barnesiella::Bacteroides 0.000 0.000
## Porphyromonas::Gelidibacter 0.000 0.000
## Bacilloplasma 0.000 0.000
## Natronincola 0.000 0.000
## Lumbricincola 0.000 0.000
## Desulfotomaculum::Clostridium 0.000 0.000
## Roseburia::Ruminococcus 0.000 0.000
## Bifidobacterium 0.000 0.000
## Streptacidiphilus 0.000 0.000
## Butyrivibrio::Pseudobutyrivibrio 0.101 0.995
## Aeromicrobium 0.000 0.000
## Proteus 0.000 0.000
## Sporobacterium::Clostridium 0.000 0.000
## Butyrivibrio::Roseburia 0.000 0.000
## Luteococcus 0.000 0.000
## Clostridium::Blautia::Desulfotomaculum 0.101 0.995
## Lachnospira::Ruminococcus 0.101 0.995
## Ornithinimicrobium 0.000 0.000
## Persicivirga 0.000 0.000
## Lachnospira::Ruminococcus::Escherichia 0.000 0.000
## Anaerophaga 0.211 -0.977
## Dysgonomonas::Flavobacterium 0.000 0.000
## Bizionia 0.101 0.995
## r2 Pr(>r)
## Tannerella 0.70 0.041 *
## Lactococcus 0.40 0.264
## Lactobacillus 0.84 0.012 *
## Lactobacillus::Lactococcus 0.35 0.327
## Parasutterella 0.17 0.674
## Helicobacter 0.57 0.050 *
## Prevotella 0.52 0.162
## Bacteroides 0.44 0.262
## Barnesiella 0.44 0.158
## Odoribacter 0.39 0.293
## Eubacterium 0.48 0.201
## Allobaculum 0.41 0.241
## Roseburia 0.15 0.807
## Clostridium 0.58 0.120
## Porphyromonas 0.63 0.040 *
## Butyrivibrio 0.30 0.432
## Ruminococcus 0.38 0.277
## Acholeplasma 0.46 0.205
## Alistipes 0.23 0.579
## Clostridium::Coprococcus 0.00 1.000
## Eubacterium (Erysipelotrichaceae)::Eubacterium 0.14 0.680
## Hydrogenoanaerobacterium 0.14 0.711
## Paraprevotella 0.56 0.081 .
## Blautia 0.03 0.939
## Adlercreutzia::Asaccharobacter 0.06 0.978
## Coprococcus 0.37 0.312
## Parabacteroides 0.50 0.097 .
## Eubacterium (Erysipelotrichaceae) 0.25 0.543
## Oscillibacter 0.34 0.356
## Acinetobacter 0.37 0.361
## Alkalitalea 0.35 0.282
## Sporobacter 0.00 1.000
## Oscillospira 0.17 0.588
## Blautia::Clostridium 0.08 0.803
## Planctomyces 0.21 0.752
## Akkermansia 0.15 0.648
## Ruminofilibacter 0.38 0.246
## Pseudoflavonifractor 0.37 0.299
## Butyricimonas 0.45 0.212
## Desulfovibrio 0.40 0.217
## Anaerotruncus 0.23 0.603
## Alistipes::Bacteroides 0.22 0.615
## Turicibacter 0.46 0.191
## Mucispirillum 0.37 0.361
## Lachnospira 0.11 0.865
## Catabacter 0.00 1.000
## Desulfosporosinus 0.00 1.000
## Bacillus 0.64 0.054 .
## Sporanaerobacter 0.22 0.629
## Streptomyces 0.10 0.868
## Butyrivibrio::Clostridium 0.39 0.242
## Candidatus Arthromitus 0.51 0.103
## Enterorhabdus 0.15 0.713
## Kurthia 0.00 1.000
## Eggerthella 0.38 0.268
## Actinocorallia 0.00 1.000
## Caldanaerocella 0.37 0.361
## Dorea 0.17 0.621
## Streptococcus 0.07 0.872
## Adlercreutzia 0.30 0.443
## Paraeggerthella 0.17 0.685
## Lachnobacterium 0.38 0.246
## TM7 (genus) 0.10 1.000
## Clostridium::Anaerostipes 0.37 0.361
## Erysipelothrix 0.25 0.475
## Formosa 0.00 1.000
## Rikenella 0.12 0.703
## Dysgonomonas 0.46 0.167
## Desulfitobacterium 0.37 0.361
## Oribacterium 0.21 0.752
## Syntrophococcus 0.00 1.000
## Ruminococcus::Clostridium 0.38 0.246
## Rhodospirillum 0.37 0.361
## Pedobacter 0.82 0.014 *
## Acetivibrio 0.46 0.242
## Clostridium::Eubacterium 0.05 0.918
## Limibacter 0.50 0.197
## Clostridium::Ruminococcus 0.00 1.000
## Coprobacillus 0.27 0.447
## Cytophaga 0.47 0.206
## Denitrobacterium 0.04 0.870
## Mycoplasma 0.57 0.118
## Roseburia::Clostridium 0.23 0.603
## Sporobacterium 0.37 0.251
## Eubacterium::Ruminococcus 0.00 1.000
## Pseudobutyrivibrio 0.10 1.000
## Robinsoniella 0.00 1.000
## Brevibacterium 0.30 0.387
## Blautia::Ruminococcus 0.26 0.526
## Pseudomonas 0.58 0.108
## Clostridium::Roseburia 0.00 1.000
## Desulfotomaculum 0.21 0.605
## Clostridium (Erysipelotrichaceae) 0.13 0.824
## Olsenella 0.06 0.897
## Azospirillum 0.00 1.000
## Oxobacter 0.00 1.000
## Asaccharobacter::Adlercreutzia 0.00 1.000
## Desulfocurvus 0.10 1.000
## Anaerostipes 0.10 1.000
## Clostridium::Ruminococcus::Desulfotomaculum::Escherichia 0.00 1.000
## Halanaerobacter 0.24 0.532
## Slackia 0.58 0.101
## Anaplasma 0.00 1.000
## Syntrophomonas 0.00 1.000
## Clostridium::Blautia 0.00 1.000
## Anaerosporobacter::Clostridium 0.00 1.000
## Clostridium::Bacteroides 0.00 1.000
## Blautia::Lactonifactor 0.00 1.000
## Parabacteroides::Bacteroides 0.21 0.752
## Enterorhabdus::Adlercreutzia 0.00 1.000
## Flavobacterium 0.43 0.233
## Clostridium::Desulfotomaculum 0.00 1.000
## Parasporobacterium 0.37 0.361
## Coprococcus::Clostridium 0.00 1.000
## Fusobacterium 0.36 0.281
## Ruminococcus::Escherichia 0.23 0.603
## Clostridium::Ruminococcus::Coprococcus 0.00 1.000
## Lactonifactor 0.00 1.000
## Haemophilus 0.06 0.832
## Sporichthya 0.26 0.513
## Cytophaga::Flavobacterium 0.38 0.268
## Clostridium::Dorea 0.00 1.000
## Ruminococcus::Blautia 0.00 1.000
## Frigoribacterium 0.21 0.752
## Paenibacillus 0.29 0.450
## Johnsonella 0.10 1.000
## Pseudoflavonifractor::Clostridium 0.00 1.000
## Adlercreutzia::Enterorhabdus::Asaccharobacter 0.00 1.000
## Nocardiopsis 0.21 0.752
## Pedobacter::Pseudomonas 0.00 1.000
## Flexibacter 0.42 0.245
## Catabacter::Ruminococcus 0.00 1.000
## Butyricimonas::Bacteroides 0.21 0.752
## Staphylococcus 0.50 0.162
## Alkalitalea::Prevotella 0.00 1.000
## Lachnospira::Anaerostipes 0.00 1.000
## Flavobacterium::Cytophaga 0.38 0.246
## Gelidibacter 0.46 0.197
## Treponema 0.21 0.752
## Pontibacter 0.00 1.000
## Desulfuromonas 0.00 1.000
## Butyricicoccus 0.00 1.000
## Clostridium::Butyrivibrio 0.23 0.603
## Coprobacillus::Clostridium 0.00 1.000
## Porphyromonas::Prevotella 0.00 1.000
## Effluviibacter 0.26 0.487
## Caminicella 0.10 1.000
## Prevotella::Bacteroides 0.00 1.000
## Pseudoalteromonas 0.00 1.000
## Butyrivibrio::Blautia 0.00 1.000
## Lachnospira::Clostridium 0.15 0.884
## Alicyclobacillus 0.00 1.000
## Lachnospira::Robinsoniella 0.00 1.000
## Subdoligranulum 0.00 1.000
## Anaerofilum 0.00 1.000
## Sporosarcina 0.00 1.000
## Alkalitalea::Sphingobacterium 0.00 1.000
## Kopriimonas 0.28 0.525
## Olivibacter 0.00 1.000
## Lachnobacterium::Coprococcus 0.00 1.000
## Fluviicola 0.56 0.122
## Peptococcus 0.00 1.000
## Ruminococcus::Roseburia 0.00 1.000
## Marinilabilia 0.00 1.000
## Catonella 0.00 1.000
## Sphingobium 0.10 1.000
## Olsenella::Streptomyces 0.00 1.000
## Terasakiella 0.00 1.000
## Roseospirillum 0.14 0.685
## Clostridium::Ruminococcus::Blautia 0.00 1.000
## Lachnospira::Roseburia 0.00 1.000
## Thalassospira 0.00 1.000
## Eubacterium::Clostridium 0.10 1.000
## Allobaculum::Eubacterium 0.00 1.000
## Blautia::Lachnospira 0.00 1.000
## Cryptobacterium 0.00 1.000
## Atopobium 0.00 1.000
## Eubacterium::Lachnospira 0.00 1.000
## Bacteroides::Alistipes 0.00 1.000
## Faecalibacterium 0.30 0.387
## Lachnospira::Coprococcus 0.00 1.000
## Microbacterium 0.00 1.000
## Zhangella::Stella 0.28 0.525
## Paludibacter 0.00 1.000
## Butyrivibrio::Ruminococcus 0.00 1.000
## Bacteroides::Lactobacillus 0.00 1.000
## Prevotella::Flavobacterium 0.15 0.884
## Ruminococcus::Dorea 0.00 1.000
## Slackia::Asaccharobacter::Adlercreutzia 0.00 1.000
## Caldicellulosiruptor 0.56 0.122
## Kordia 0.00 1.000
## Bacteroides::Lactobacillus::Lachnospira 0.23 0.603
## Sphingobacterium 0.00 1.000
## Anaeroplasma 0.23 0.603
## Atopostipes 0.00 1.000
## Enterococcus 0.00 1.000
## Insolitispirillum 0.00 1.000
## Clostridium::Acetivibrio 0.00 1.000
## Plantactinospora 0.00 1.000
## Stappia 0.00 1.000
## Anaplasma::Clostridium 0.00 1.000
## Clostridium (Erysipelotrichaceae)::Clostridium 0.00 1.000
## Algoriphagus 0.00 1.000
## Dorea::Ruminococcus 0.00 1.000
## Roseburia::Lachnospira 0.00 1.000
## Rhizobium 0.00 1.000
## Anaerofustis 0.00 1.000
## Echinicola 0.00 1.000
## Anaerostipes::Clostridium 0.00 1.000
## Lachnospira::Blautia 0.00 1.000
## Aestuariimicrobium 0.00 1.000
## Gelidibacter::Flavobacterium 0.00 1.000
## Helicobacter::Flexispira 0.10 1.000
## Eubacterium (Erysipelotrichaceae)::Paenibacillus::Eubacterium 0.00 1.000
## Pelospora 0.00 1.000
## Stella 0.00 1.000
## Methylobacterium 0.00 1.000
## Rickettsia 0.23 0.603
## Porphyromonas::Flavobacterium 0.00 1.000
## Adlercreutzia::Enterorhabdus 0.00 1.000
## Ruminococcus::Escherichia::Parabacteroides 0.00 1.000
## Limibacter::Bacteroides 0.00 1.000
## Paraprevotella::Prevotella 0.00 1.000
## Janthinobacterium::Zoogloea::Duganella 0.00 1.000
## Flavobacterium::Gelidibacter 0.00 1.000
## Barnesiella::Bacteroides 0.00 1.000
## Porphyromonas::Gelidibacter 0.00 1.000
## Bacilloplasma 0.00 1.000
## Natronincola 0.00 1.000
## Lumbricincola 0.00 1.000
## Desulfotomaculum::Clostridium 0.00 1.000
## Roseburia::Ruminococcus 0.00 1.000
## Bifidobacterium 0.00 1.000
## Streptacidiphilus 0.00 1.000
## Butyrivibrio::Pseudobutyrivibrio 0.10 1.000
## Aeromicrobium 0.00 1.000
## Proteus 0.00 1.000
## Sporobacterium::Clostridium 0.00 1.000
## Butyrivibrio::Roseburia 0.00 1.000
## Luteococcus 0.00 1.000
## Clostridium::Blautia::Desulfotomaculum 0.10 1.000
## Lachnospira::Ruminococcus 0.10 1.000
## Ornithinimicrobium 0.00 1.000
## Persicivirga 0.00 1.000
## Lachnospira::Ruminococcus::Escherichia 0.00 1.000
## Anaerophaga 0.56 0.122
## Dysgonomonas::Flavobacterium 0.00 1.000
## Bizionia 0.10 1.000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## Permutation: free
## Number of permutations: 999Los géneros con un \(p-\)valor menor a \(0.05\) se muestran con el siguiente gráfico.
# Correr todo junto
# Definir los parámetros del plot
par(mar = c(5, 5, 1, 2) + 0.1)
# Inicializar el plot
plot(PCoA$points[ ,1], PCoA$points[ ,2], ylim = c(-0.5, 0.5),
xlab = paste("PCoA 1 (", explainedvar1, "%)", sep = ""),
ylab = paste("PCoA 2 (", explainedvar2, "%)", sep = ""),
pch = 5, cex = 1.0, type = "n", cex.lab = 1.0, cex.axis = 1.2, axes = FALSE)
# Añadir ejes
axis(side = 1, labels = T, lwd.ticks = 2, cex.axis = 1.2, las = 1)
axis(side = 2, labels = T, lwd.ticks = 2, cex.axis = 1.2, las = 1)
abline(h = 0, v = 0, lty = 3)
box(lwd = 2)
# Añadir puntos y etiquetas
points(PCoA$points[ ,1], PCoA$points[ ,2],
pch = 19, cex = 3, bg = "blue", col = "blue")
text(PCoA$points[ ,1], PCoA$points[ ,2],
labels = row.names(PCoA$points))
fecalREL <- abund_table
for(i in 1:nrow(abund_table)){
fecalREL[i, ] = abund_table[i, ] / sum(abund_table[i, ])
}
#install.packages("BiodiversityR")
library("pbkrtest")
library("BiodiversityR")
# Calcular y añadir el score de las especies
PCoA <- add.spec.scores(PCoA,fecalREL,method = "pcoa.scores",Rscale=TRUE,scaling=1, multi=1)
text(PCoA$cproj[ ,1], PCoA $cproj[ ,2],
labels = row.names(PCoA$cproj),cex=0.5, col = "blue")
genus_corr <- add.spec.scores(PCoA, fecalREL, method = "cor.scores")$cproj
corrcut <- 0.7 # definir el cutoff
import_genus <- genus_corr[abs(genus_corr[, 1]) >= corrcut | abs(genus_corr[, 2]) >= corrcut, ]
# Código para los generos con un p-valor<0.05
fit <- envfit(PCoA, fecalREL, perm = 999)
plot(fit, p.max = 0.05, cex=0.5, col = "red")
Ordination plot para los dos ejes de PCoA de influencia de genero con un \(p-\)valor \(<0.05\)
Escalamiento multidimensional no métrico (NMDS)
El NMDS es la alternativa no métrica al análisis PCoA. El NMDS ha sido reconocido como un buen método de ordenación porque utiliza formas ecológicamente significativas para medir las disimilitudes de las comunidades y cualquier medida de distancia (disimilitudes) entre las muestras como entrada. Por ello, es el método recomendado en la ordenación de comunidades. El objetivo principal del análisis NMDS es proyectar la posición relativa de los puntos de la muestra en un espacio de ordenación de baja dimensión (dos o tres ejes).
La función metaMDS() del paquete vegan
realiza el análisis NMDS. Para simplificar, el algoritmo del análisis
NMDS se resume como sigue: en primer lugar, utiliza la función
vegdist() para obtener las medidas de disimilitud
adecuadas; a continuación, ejecuta NMDS varias veces con configuraciones
iniciales aleatorias y compara los resultados mediante la función
procrustes(), y se detiene tras encontrar dos veces una
solución mínima similar. Por último, escala y rota la solución, y añade
las puntuaciones de las especies (o OTUs/taxa) a la configuración como
promedios ponderados utilizando la función wascores(). Una
vez finalizado el algoritmo, la solución final se rota utilizando PCA
para facilitar su interpretación.
En este caso, utilizamos el paquete vegan y los mismos
datos fecales de ratones Vdr-/- para ilustrar análisis NMDS. En primer
lugar, utilizamos la función vegdist() para obtener las
medidas de disimilitud de Bray-Curtis (el método por defecto)
utilizando la configuración por defecto de la función
metaMDS(). Esta función transforma automáticamente los
datos y comprueba la solidez de la solución. La doble estandarización de
Wisconsin y la transformación sqrt se
utilizaron conjuntamente en la llamada a la función
metaMDS(). La llamada a la función produce un valor de
tensión del 7,57%.
bc_nmds <- metaMDS(abund_table, dist = "bray")## Square root transformation
## Wisconsin double standardization
## Run 0 stress 0.07574
## Run 1 stress 0.07574
## ... New best solution
## ... Procrustes: rmse 6.721e-07 max resid 1.058e-06
## ... Similar to previous best
## Run 2 stress 0.1192
## Run 3 stress 0.1284
## Run 4 stress 0.07574
## ... Procrustes: rmse 3.801e-07 max resid 6.028e-07
## ... Similar to previous best
## Run 5 stress 0.1758
## Run 6 stress 0.1162
## Run 7 stress 0.1192
## Run 8 stress 0.07574
## ... Procrustes: rmse 1.075e-06 max resid 1.813e-06
## ... Similar to previous best
## Run 9 stress 0.1284
## Run 10 stress 0.1284
## Run 11 stress 0.1899
## Run 12 stress 0.1393
## Run 13 stress 0.1562
## Run 14 stress 0.1801
## Run 15 stress 0.229
## Run 16 stress 0.1562
## Run 17 stress 0.1192
## Run 18 stress 0.1284
## Run 19 stress 0.1192
## Run 20 stress 0.1162
## *** Solution reachedbc_nmds##
## Call:
## metaMDS(comm = abund_table, distance = "bray")
##
## global Multidimensional Scaling using monoMDS
##
## Data: wisconsin(sqrt(abund_table))
## Distance: bray
##
## Dimensions: 2
## Stress: 0.07574
## Stress type 1, weak ties
## Two convergent solutions found after 20 tries
## Scaling: centring, PC rotation, halfchange scaling
## Species: expanded scores based on 'wisconsin(sqrt(abund_table))'En segundo lugar, utilizamos la función ordiplot() para
dibujar los resultados del NMDS. La configuración por defecto sólo añade
puntos a la figura, utilizamos type = 't' o
type = 'text' para añadir etiquetas de texto en esta
figura.
ordiplot (bc_nmds, type = 't')
Ordiplot de dos ejes de NMDS con etiquetas
Para trazar las puntuaciones del sitio (muestra) como texto:
ordiplot(bc_nmds, display = "sites", type = "text")
Ordiplot de dos ejes de NMDS con etiquetas de las muestras
Finalmente, la función stressplot() se utiliza para
dibujar el gráfico de estrés de Shepards. La función
stressplot() genera dos figuras: una traza las distancias
de ordenación frente a la disimilitud observada (las disimilitudes de la
comunidad elegida), junto con una línea de paso monótona para mostrar el
ajuste; otra traza la bondad del ajuste para evaluar la bondad de la
ordenación de NMDS2 frente a NMDS1 de muestras concretas. Primero,
dividimos la ventana de trazado en dos paneles. La función
plot() dibuja el diagrama de ordenación NMDS con los sitios
(muestras). La función points() añade los puntos con un
tamaño que refleja la bondad del ajuste (mayor = peor ajuste).
par (mfrow = c(1,2))
stressplot (bc_nmds)
plot (bc_nmds, display = 'sites', type = 't', main = 'Goodness of fit')
points (bc_nmds, display = 'sites', cex = goodness (bc_nmds)*300)
Gráfico del strés de Shepards y de la bondad de ajuste
El stressplot muestra la relación entre las distancias reales entre las muestras en solución de ordenación de dimensión \(m\), y sus disimilitudes composicionales particulares expresadas por la medida de disimilitud de Bray-Curtis. Existen dos estadísticas de bondad de ajuste similares a la correlación; la correlación basada en la tensión \(R^2 = 1 - S^2\) (ajuste no métrico \(= 0.994\)) y la correlación entre los valores ajustados y las distancias de ordenación, o entre la línea de paso y los puntos: “\(R^2\) basado en el ajuste” (ajuste lineal \(= 0,956\)).
Análisis de correspondencia (AC)
CA es otro método de ordenación sin restricciones. Utiliza las distancias chi-cuadrado entre las muestras en el espacio multidimensional de todos los ejes de ordenación y da un alto peso a las especies raras (por ejemplo, especies de baja ocurrencia con muchos ceros). Al igual que PCA y PCoA, CA utiliza los valores propios para medir la importancia de los ejes ortogonales devueltos.
En el diagrama de ordenación del PCA, los taxones son vectores y las
muestras son puntos, mientras que en el CA, los taxones y las muestras
están representados por puntos. Al igual que PCA, CA tiene dos tipos de
escalas: Escalas 1 y 2. En el espacio de ordenación reducido, las
distancias entre las muestras (Escala 1) se aproximan a su
distancia chi-cuadrado. Por ejemplo, cualquier muestra cercana
al punto que representa un taxón probablemente contiene una alta
contribución a ese taxón. Las distancias entre taxones (Escala
2) también se aproximan a sus distancias chi-cuadrado. Por
ejemplo, cualquier taxón cercano al punto que representa una muestra
probablemente tenga una mayor frecuencia en esa muestra. La función
cca() del paquete vegan puede utilizarse para
realizar un análisis de correspondencia sin restricciones. Utilizamos la
función cca() para realizar el CCA y la función
evplot() para seleccionar ejes de ordenación importantes
basados en el criterio de Kaiser-Guttman o modelo de
broken-stick. La función evplot(), escrita por
Borcard et al. (2011), se utiliza aquí para trazar los valores propios y
los porcentajes de variación de un objeto de ordenación. Cuando se
utiliza la función cca() (cca = análisis canónico
compoente) para realizar CA sin restricciones, no se especifica la
matriz ambiental ni la información de agrupación.
fecal_genus_cca=cca(abund_table)
fecal_genus_cca## Call: cca(X = abund_table)
##
## Inertia Rank
## Total 0.502
## Unconstrained 0.502 7
## Inertia is scaled Chi-square
## 122 species (variables) deleted due to missingness
##
## Eigenvalues for unconstrained axes:
## CA1 CA2 CA3 CA4 CA5 CA6 CA7
## 0.2036 0.1256 0.0747 0.0374 0.0322 0.0200 0.0090La heterogeneidad total de los datos (inercia) es de 0.502, y el primer eje capta 40.5638% de la variación total en la composición del género (0.2036 / 0.502 = 0.4056, donde 0.2036 es el valor propio del primer eje CA1, y 0.502 es la heterogeneidad total de los datos).
Los siguientes códigos de R se utilizan para trazar la ordenación y mostrar los nombres de las muestras en la figura:
plot(fecal_genus_cca, display="sites") 
Plot de ordenación con muestras.
Si no quiere que el gráfico esté sobrecargado, utilice los siguientes códigos de R para sólo mostrar puntos en lugar de nombres de muestras en la figura:
plot(fecal_genus_cca, display="sites", type="p")
Plot de ordenación sin etiquetas de las muestras
El siguiente diagrama de ordenación revela el patrón de muestras y géneros en diagrama de ordenación:
ordiplot(fecal_genus_cca)
Plot de ordenación, muestra los patrones de las muestras y generos en el diagrama de ordenación. Los círculos representan las muestras y los signos + representan los géneros.
El siguiente análisis posterior a la CA utilizando la función
evplot() es para ilustrar cómo decidir qué eje de CA debe
utilizarse para la interpretación de los resultados. La sintaxis de la
función función evplot() es evplot(ev), donde
ev es un vector de valores propios. En primer lugar, ejecutamos la
función.
evplot <- function(ev)
{# Broken stick model (MacArthur 1957)
n <- length(ev)
bsm <- data.frame(j=seq(1:n), p=0)
bsm$p[1] <- 1/n
for (i in 2:n) bsm$p[i] <- bsm$p[i-1] + (1/(n + 1 - i))
bsm$p <- 100*bsm$p/n
# Graficar eigenvalores y porcentaje de varianza en cada eje
op <- par(mfrow=c(2,1))
barplot(ev, main="Eigenvalues", col="bisque", las=2)
abline(h=mean(ev), col="red")
legend("topright", "Average eigenvalue", lwd=1, col=2, bty="n")
barplot(t(cbind(100*ev/sum(ev), bsm$p[n:1])), beside=TRUE,
main="% variation", col=c("bisque",2), las=2)
legend("topright", c("% eigenvalue", "Broken stick model"),
pch=15, col=c("bisque",2), bty="n")
par(op)
}
# Plot de eigenvalores y porcentaje de varianza en cada eje
ev <- fecal_genus_cca$CA$eig
windows(title="CA eigenvalues")
evplot(ev)
Eigenvalores y porcentajes de variación del criterio Keiser-Guttman y el modelo broken-stick
Como se presenta en la sección de PCoA, para evaluar si los primeros ejes de CA del análisis capturan una cantidad desproporcionadamente grande de la variación total explicada, dos vías adicionales de uso de gráficos son el criterio de Keiser-Guttman y **criterio de broken-stick. Según el criterio de Kaiser-Guttman, los valores propios asociados a los primeros ejes deben ser mayores que la media de todos los valores propios; con el modelo de broken-stick, los valores propios asociados a los primeros ejes se comparan con las expectativas. Tanto el criterio de Keiser-Guttman* como el modelo de broken-stick muestran que los tres primeros ejes son importantes.
Ordenaciones restringidas
Como contrapartida a los tres métodos básicos de ordenación, el
paquete vegan cuenta con tres versiones de ordenación
restringida:
- Análisis de redundancia (RDA, relacionado con el análisis de componentes principales)
- Análisis restringido de coordenadas principales (CAP, relacionado con el escalado métrico)
- Análisis restringido de correspondencia (CCA, relacionado con el análisis de correspondencia).
Análisis de redundancia (RDA)
El RDA en la función rda() se basa en las distancias
euclidianas y combina la regresión múltiple con el PCA. Con RDA, cada
eje canónico es una combinación lineal de todas las variables
explicativas. El número de ejes canónicos corresponde al número de
variables explicativas, o más exactamente al número de grados de
libertades (Borcard et al. 2011).
Puede utilizar dos sintaxis diferentes para calcular una RDA mediante
la función rda() de vegan: sintaxis matricial
y sintaxis de fórmula. La sintaxis matricial es la más sencilla: basta
con enumerar los nombres de los objetos separados por comas como se
indica a continuación. RDA = rda(Y, X, W), donde Y es la
matriz de respuesta (composición de taxones), X es la matriz explicativa
(factores ambientales) y W es la matriz opcional de covariables.
La sintaxis de fórmula esta dada como:
RDA = rda(Y ~ var1 + factor A + var2 * var3 + condition(var4), data = XW)
donde Y es la matriz de respuesta (composición de
taxones); y las variables del lado derecho de la fórmula son las
variables explicatorias o restringidas. El objeto XW es un
dataframe con las variables explicatorias o covariables.
La prueba de hipótesis no es el caso en PCA. Sin embargo, con dos
conjuntos de datos \(Y\) y \(X\), en RDA podemos probar una hipótesis
nula de ausencia de relación lineal entre entre ellos. Utilizamos el
conjunto de datos de fumadores para ilustrar el uso de RDA. Primero,
cargue los datos del paquete GUniFrac y utilice la función
select() del paquete dplyr para crear un
subconjunto de variables explicativas.
library("GUniFrac")
data(throat.otu.tab)
data(throat.meta)
library("dplyr")##
## Attaching package: 'dplyr'## The following object is masked from 'package:kableExtra':
##
## group_rows## The following objects are masked from 'package:igraph':
##
## as_data_frame, groups, union## The following objects are masked from 'package:stats':
##
## filter, lag## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, unionthroat_meta <- select(throat.meta, SmokingStatus, Age, Sex, PackYears)A efectos de RDA, transformamos los datos de los taxones utilizando la transformación de Hellinger:
library("vegan") # Solo si no se ha cargado antes
abund_hell <- decostand (throat.otu.tab, 'hell') La RDA se calcula con la función rda() si se suministra
la matriz de variables ambientales (si no, se calculará el ACP)
rda_hell<- rda(abund_hell ~ ., throat_meta)
head(summary (rda_hell))##
## Call:
## rda(formula = abund_hell ~ SmokingStatus + Age + Sex + PackYears, data = throat_meta)
##
## Partitioning of variance:
## Inertia Proportion
## Total 0.4627 1.000
## Constrained 0.0482 0.104
## Unconstrained 0.4145 0.896
##
## Eigenvalues, and their contribution to the variance
##
## Importance of components:
## RDA1 RDA2 RDA3 RDA4 PC1 PC2 PC3
## Eigenvalue 0.0239 0.0133 0.00597 0.00504 0.0683 0.0595 0.0355
## Proportion Explained 0.0517 0.0287 0.01290 0.01089 0.1476 0.1285 0.0768
## Cumulative Proportion 0.0517 0.0805 0.09335 0.10424 0.2518 0.3803 0.4572
## PC4 PC5 PC6 PC7 PC8 PC9 PC10
## Eigenvalue 0.024 0.0174 0.0156 0.0149 0.0122 0.0103 0.00965
## Proportion Explained 0.052 0.0376 0.0337 0.0322 0.0263 0.0223 0.02086
## Cumulative Proportion 0.509 0.5467 0.5804 0.6126 0.6390 0.6613 0.68212
## PC11 PC12 PC13 PC14 PC15 PC16 PC17
## Eigenvalue 0.00924 0.00876 0.00775 0.00671 0.00648 0.00626 0.00597
## Proportion Explained 0.01998 0.01893 0.01674 0.01450 0.01399 0.01352 0.01291
## Cumulative Proportion 0.70210 0.72103 0.73777 0.75226 0.76626 0.77978 0.79269
## PC18 PC19 PC20 PC21 PC22 PC23 PC24
## Eigenvalue 0.00549 0.00529 0.0050 0.00459 0.00435 0.00400 0.00372
## Proportion Explained 0.01187 0.01144 0.0108 0.00993 0.00939 0.00865 0.00804
## Cumulative Proportion 0.80455 0.81599 0.8268 0.83673 0.84612 0.85477 0.86281
## PC25 PC26 PC27 PC28 PC29 PC30 PC31
## Eigenvalue 0.00361 0.00357 0.00338 0.00323 0.00314 0.00307 0.00296
## Proportion Explained 0.00780 0.00771 0.00730 0.00698 0.00678 0.00663 0.00639
## Cumulative Proportion 0.87062 0.87832 0.88562 0.89260 0.89938 0.90601 0.91240
## PC32 PC33 PC34 PC35 PC36 PC37 PC38
## Eigenvalue 0.00285 0.00268 0.00249 0.00232 0.00225 0.00218 0.00213
## Proportion Explained 0.00617 0.00579 0.00537 0.00502 0.00486 0.00471 0.00460
## Cumulative Proportion 0.91857 0.92435 0.92973 0.93475 0.93961 0.94431 0.94891
## PC39 PC40 PC41 PC42 PC43 PC44 PC45
## Eigenvalue 0.00198 0.00186 0.00176 0.00168 0.00163 0.00155 0.00151
## Proportion Explained 0.00427 0.00403 0.00380 0.00364 0.00353 0.00335 0.00327
## Cumulative Proportion 0.95319 0.95721 0.96101 0.96465 0.96818 0.97153 0.97480
## PC46 PC47 PC48 PC49 PC50 PC51 PC52
## Eigenvalue 0.00145 0.00139 0.00135 0.00123 0.00121 0.00117 0.00103
## Proportion Explained 0.00313 0.00301 0.00292 0.00266 0.00261 0.00253 0.00222
## Cumulative Proportion 0.97793 0.98094 0.98386 0.98652 0.98913 0.99165 0.99387
## PC53 PC54 PC55
## Eigenvalue 0.00100 0.000943 0.000892
## Proportion Explained 0.00216 0.002037 0.001927
## Cumulative Proportion 0.99604 0.998073 1.000000
##
## Accumulated constrained eigenvalues
## Importance of components:
## RDA1 RDA2 RDA3 RDA4
## Eigenvalue 0.0239 0.0133 0.00597 0.00504
## Proportion Explained 0.4961 0.2757 0.12374 0.10446
## Cumulative Proportion 0.4961 0.7718 0.89554 1.00000
##
## Scaling 2 for species and site scores
## * Species are scaled proportional to eigenvalues
## * Sites are unscaled: weighted dispersion equal on all dimensions
## * General scaling constant of scores: 2.286
##
##
## Species scores
##
## RDA1 RDA2 RDA3 RDA4 PC1 PC2
## 4695 0.001921 -0.011290 -0.007088 -0.002362 -0.000786 -0.000549
## 2983 0.002385 -0.001850 -0.004702 0.004402 0.000896 -0.003710
## 2554 0.001181 -0.005402 -0.004723 -0.001745 -0.005643 -0.004346
## 3315 -0.000150 -0.000566 -0.001440 -0.000609 -0.000684 -0.007786
## 879 -0.003149 -0.002243 0.000376 -0.002867 -0.000290 0.003117
## 1313 0.000261 0.003148 0.003128 0.000805 0.000496 -0.009699
## ....
##
##
## Site scores (weighted sums of species scores)
##
## RDA1 RDA2 RDA3 RDA4 PC1 PC2
## ESC_1.1_OPL -0.444 -0.4795 -0.1746 -0.3520 0.0816 -0.1354
## ESC_1.3_OPL 0.601 -0.2263 -0.1055 -0.2759 0.1837 -0.0773
## ESC_1.4_OPL 0.657 0.1918 -0.1318 0.2111 0.3182 0.1287
## ESC_1.5_OPL 0.821 0.0243 -0.0463 0.5609 0.4416 0.3979
## ESC_1.6_OPL 0.980 -0.3117 0.4864 -0.0419 0.3708 -0.0187
## ESC_1.10_OPL 0.117 0.5722 0.8475 -0.2146 -0.3048 0.3063
## ....
##
##
## Site constraints (linear combinations of constraining variables)
##
## RDA1 RDA2 RDA3 RDA4 PC1 PC2
## ESC_1.1_OPL -0.3062 -0.22035 0.0309 -0.3282 0.0816 -0.1354
## ESC_1.3_OPL 0.3242 0.05576 -0.1321 0.1074 0.1837 -0.0773
## ESC_1.4_OPL 0.3073 0.06338 -0.1988 0.1493 0.3182 0.1287
## ESC_1.5_OPL 0.0932 -0.50693 -0.3499 0.0537 0.4416 0.3979
## ESC_1.6_OPL 0.3660 0.00406 0.1952 0.0548 0.3708 -0.0187
## ESC_1.10_OPL 0.4675 -0.10112 0.8899 -0.1050 -0.3048 0.3063
## ....
##
##
## Biplot scores for constraining variables
##
## RDA1 RDA2 RDA3 RDA4 PC1 PC2
## SmokingStatusSmoker 0.906 -0.3399 -0.135 0.2110 0 0
## Age 0.228 -0.0442 0.748 0.6217 0 0
## SexMale 0.496 0.8360 0.103 0.2110 0 0
## PackYears 0.658 -0.1515 0.738 0.0222 0 0
##
##
## Centroids for factor constraints
##
## RDA1 RDA2 RDA3 RDA4 PC1 PC2
## SmokingStatusNonSmoker -0.250 0.0938 0.0374 -0.0582 0 0
## SmokingStatusSmoker 0.286 -0.1072 -0.0427 0.0665 0 0
## SexFemale -0.199 -0.3362 -0.0413 -0.0848 0 0
## SexMale 0.107 0.1810 0.0222 0.0457 0 0Los cuatro ejes de ordenación restringidos (denominados RDA1 a RDA4) están relacionados con las 4 variables ambientales (SmokingStatus, Age, Sex, PackYears). Los ejes no restringidos se denominan ejes PC.
Las varianzas totales se dividen en varianza restringida y varianza no restringida. La varianza restringida se explica por los ejes restringidos (es decir, las variables ambientales); mientras que la varianza no restringida se explica por los ejes no restringidos (es decir, la varianza no explicada por factores ambientales). La varianza restringida es la cantidad de varianza que la matriz \(Y\) es explicada por las variables explicativas. Es un equivalente a un R2 sesgado y no ajustado en la regresión múltiple. La tabla de partición de la varianza muestra que el 4.8% (0.0482) de la proporción de la varianza total es explicada por estos 4 factores ambientales.
La función coef() recupera los coeficientes canónicos
(el equivalente a coeficientes de regresión) para cada variable
explicativa en cada eje canónico.
coef(rda_hell)## RDA1 RDA2 RDA3 RDA4
## SmokingStatusSmoker 0.177040 -0.120402 -0.187134 0.15184
## Age -0.004452 -0.002161 0.002997 0.01757
## SexMale 0.092005 0.267873 -0.027579 0.01856
## PackYears 0.005257 -0.000520 0.011633 -0.01596Podemos utilizar la función RsquareAdj() para extraer el
valor de \(R^2\) y el ajustado \(R^2\) de los resultados de la
ordenación.
RsquareAdj(rda_hell)## $r.squared
## [1] 0.1042
##
## $adj.r.squared
## [1] 0.03909La función devuelve dos \(R^2\): uno es el \(R^2\) ordinario (r.al cuadrado), otro es el \(R^2_{adj}\) ajustado. Tenga en cuenta que \(R^2_{adj}\) es siempre menor que el \(R^2\), y la diferencia aumenta con el aumento del número de variables explicativas. El \(R^2_{adj}\) puede ser negativo, lo que significa que las variables explicativas explican menos variación que el mismo número de variables generadas al azar.
Aplicamos el criterio de Kaiser-Guttman a los ejes residuales.
rda_hell$CA$eig[rda_hell$CA$eig > mean(rda_hell$CA$eig)]## PC1 PC2 PC3 PC4 PC5 PC6 PC7 PC8
## 0.068299 0.059456 0.035542 0.024040 0.017406 0.015601 0.014898 0.012182
## PC9 PC10 PC11 PC12 PC13
## 0.010314 0.009654 0.009245 0.008760 0.007745A continuación, trazamos los resultados de la RDA. En función de la ordenación que le interese, puede elegir trazar el escalado tipo 1 o el escalado tipo 2. Para la ordenación de las muestras, el escalamiento 1 es la opción más adecuada. Para la ordenación de los taxones, entonces, el escalamiento 2 es el más apropiado.
El escalamiento de tipo 1 hace hincapié en las relaciones entre las muestras. Las características esenciales (Borcard et al. 2011; Legendre y Legendre 2012) son que las muestras actúan como los centroides de las variables de respuesta (columnas) y las distancias entre los puntos de las muestras indican sus distancias \(\chi^2\). La interpretación se basa en que:
- Los puntos de muestra que están cerca unos de otros probablemente sean relativamente similares en cuanto a sus frecuencias relativas.
- Las muestras cercanas a los centroides que representan estados de variables categóricas o cualitativas tienen más probabilidades de poseer el estado para esa especie/taxón.
- Una proyección en ángulo recto de un punto de muestra sobre un vector que representa una variable cuantitativa cuantitativa se aproxima al valor de la variable realizada para esa muestra.
El escalamiento de tipo 2 hace hincapié en las relaciones entre las variables de respuesta. Las características esenciales son que las variables de respuesta (columnas) actúan como los centroides de las muestras y las distancias entre los puntos de las variables de respuesta indican sus distancias \(\chi^2\). La interpretación se basa en que:
- Los puntos de especies/taxones que están cerca unos de otros entre sí es probable que tengan frecuencias relativas similares a lo largo de las muestras.
- Cuanto más cerca una variable de respuesta (una especie/taxón) al centroide que representa un estado de una categórica, es más probable que esa variable de respuesta tenga valores más altos en ese estado.
- Una proyección en ángulo recto de un punto que representa una respuesta (una especie/taxón) sobre una flecha que representa una variable explicativa indica la posición del valor máximo (el óptimo) de la variable de respuesta a lo largo de esa variable explicativa.
En la ordenación restringida, se dispone de una fuente de datos adicional para el trazado: las variables explicativas (ambientales). Así, se dispone de tres datos diferentes para muestras, variables de respuesta y variables explicativas. Si se muestran todos los datos de las tres fuentes, se tiene un “triplot”. Si muestra los datos de dos fuentes, tiene biplot. Los siguientes códigos de R generan un triplot de escala de tipo 2:
plot(rda_hell, display=c("sp", "lc", "cn"),main="Triplot RDA - scaling 2")
taxa_scores <- scores(rda_hell, choices=c(1,2), display="sp")
arrows(0, 0, taxa_scores[,1], taxa_scores[,2], length=0, lty=1, col="red")
Triplot RDA con escalado tipo 2
En los códigos anteriores, “sp” representa las especies
(display = "sp"), y “cn” las restricciones o las variables
explicativas (display = "cn"). Hay dos puntuaciones
muestrales en el paquete vegan: sumas ponderadas de
especies/taxa (display = "wa"), y puntuaciones muestrales
ajustadas o “LC scores” (display = "lc"), es decir,
combinaciones lineales de variables explicativas. Se puede elegir una de
ellas para mostrarla en un gráfico. La segunda línea de código de R se
utilizan para añadir flechas para mostrar especies/taxa. El argumento de
la función scores() "choices=c(1,2)" es los
ejes a trazar. Los valores por defecto trazan los ejes 1 y 2. Por lo
tanto, podemos omitir la opción "choices=c(1,2)". En el
gráfico \(\ref{fig:triplot_RDA}\), los
triángulos huecos verdes representan muestras, las cruces azules
representan los estados de una variable explicativa categórica (por
ejemplo, hombre, mujer o fumador, no fumador), y las flechas azules para
las variables explicativas cuantitativas (en este caso, PackYears y
Edad) con puntas de flecha que indican su dirección de aumento, y las
especies/taxas se muestran en rojo.
Para comprobar la importancia de la variación de los datos de la
comunidad de fumadores explicada por las variables explicativas, podemos
realizar una prueba de permutación de Monte Carlo mediante la función
función anova.cca() del paquete vegan. Esta
función puede probar la significancia del modelo global (por defecto),
todos los ejes (by = "axis"), las variables explicativas
individuales (by = "terms"), el primer eje restringido
(first = TRUE), o la variación explicada por variables
explicativas individuales después de eliminar la variación de todas las
demás variables en el modelo (by = "margin"). Ilustramos su
capacidad probando el modelo global, cada eje y cada variable
explicativa, respectivamente.
Para asegurarnos de obtener el mismo resultado de la prueba de permutación cada vez que se ejecute, establecemos primero la misma semilla. Los siguientes códigos R prueban la significación del modelo global. El argumento “paso” especifica el número mínimo de permutaciones.
set.seed (123)
anova(rda_hell, step=1000) %>%
kbl(booktabs = TRUE, align = "c",
caption = "ANOVA global")%>%
kable_styling(position = "center",
latex_options = c("hold_position", "striped"),
font_size = 12)| Df | Variance | F | Pr(>F) | |
|---|---|---|---|---|
| Model | 4 | 0.0482 | 1.6 | 0.004 |
| Residual | 55 | 0.4145 | NA | NA |
Podemos observar que el test global es estadísticamente significante. Ahora, vamos a ver el test en cada eje.
set.seed (123)
anova(rda_hell, by="axis", step=1000)%>%
kbl(booktabs = TRUE, align = "c",
caption = "ANOVA test para todos los ejes")%>%
kable_styling(position = "center",
latex_options = c("hold_position", "striped"),
font_size = 12)| Df | Variance | F | Pr(>F) | |
|---|---|---|---|---|
| RDA1 | 1 | 0.0239 | 3.1754 | 0.013 |
| RDA2 | 1 | 0.0133 | 1.7644 | 0.256 |
| RDA3 | 1 | 0.0060 | 0.7920 | 0.942 |
| RDA4 | 1 | 0.0050 | 0.6686 | 0.876 |
| Residual | 55 | 0.4145 | NA | NA |
Podemos observar que el primer y segundo eje son significantes. El siguiente código es el test de significancia para cada variable explicativa.
set.seed (123)
anova(rda_hell, by="terms", step=1000)%>%
kbl(booktabs = TRUE, align = "c",
caption = "ANOVA test para variables explicativas")%>%
kable_styling(position = "center",
latex_options = c("hold_position", "striped"),
font_size = 12)| Df | Variance | F | Pr(>F) | |
|---|---|---|---|---|
| SmokingStatus | 1 | 0.0215 | 2.8571 | 0.002 |
| Age | 1 | 0.0057 | 0.7512 | 0.773 |
| Sex | 1 | 0.0144 | 1.9171 | 0.019 |
| PackYears | 1 | 0.0066 | 0.8751 | 0.588 |
| Residual | 55 | 0.4145 | NA | NA |
Aquí podemos observar que tanto SmokingStatus como
Sex son significativas.
Por último, utilizamos la selección directa para reducir el número de variables explicativas que entran en el análisis, mientras se optimiza la variación explicada por ellas.
Actualmente, hay tres funciones disponibles en RDA para la selección
hacia adelante: ordistep(), ordiR2step(), y
forward.sel().Las dos primeras funciones son del paquete
vegan. La tercera, está disponible en el paquete
adespatial. Estas funciones utilizan diferentes criterios
para la selección de variables. La función ordistep()
utiliza el criterio AIC y los \(p-\)valores de la prueba de permutación de
Monte Carlo para la comparación de variables. La función
ordiR2step() utiliza \(R_{adj}^2\). La función
forward.sel() utiliza el nivel de significación
preseleccionado de \(\alpha\) como
criterio de selección de variables.
Aunque la función forward.sel() tiene una lógica
diferente para establecer los argumentos comparada con las otras dos
funciones, básicamente los resultados devueltos son los mismos que
ordiR2step(). Por lo tanto, sólo ilustramos las dos
primeras funciones. El procedimiento ordistep() es
aplicable con las funciones rda(), cca() o
cmdscale(). El procedimiento ordiR2step() sólo
puede aplicarse a rda() y capscale(), pero no
para cca() porque cca() no devuelve \(R_{adj}^{2}\).
Los siguientes códigos de R utilizan la función
ordistep() para ejecutar la selección hacia adelante.
También se dispone de opciones para la selección por pasos y hacia
atrás. Esta función permite el uso de factores.
step_forward <- ordistep(rda(abund_hell ~ 1, data=throat_meta),
scope=formula(rda_hell), direction="forward", pstep=1000)##
## Start: abund_hell ~ 1
##
## Df AIC F Pr(>F)
## + SmokingStatus 1 -46.1 2.83 0.005 **
## + Sex 1 -45.3 2.01 0.015 *
## + PackYears 1 -45.1 1.80 0.050 *
## + Age 1 -44.1 0.83 0.570
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Step: abund_hell ~ SmokingStatus
##
## Df AIC F Pr(>F)
## + Sex 1 -46.0 1.87 0.04 *
## + PackYears 1 -45.0 0.87 0.63
## + Age 1 -44.9 0.74 0.78
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Step: abund_hell ~ SmokingStatus + Sex
##
## Df AIC F Pr(>F)
## + PackYears 1 -45.0 0.85 0.66
## + Age 1 -44.9 0.81 0.66Este procedimiento elige las variables SmokingStatus y
Sex. La siguiente selección utiliza la función
ordiR2step(). La configuración por defecto de este
procedimiento para incluir una nueva variable se basa en \(R^2_{adj}\) y su comparación con la del
modelo global (con todas las variables). La selección se detendrá si la
nueva variable no es significativa o el \(R_{adj}^2\) del modelo que incluye esta
nueva variable superaría el \(R_{adj}^{2}\) del modelo global.
step_forward <- ordiR2step(rda(abund_hell ~ 1, data=throat_meta),
scope=formula(rda_hell), direction="forward", pstep=1000)## Step: R2.adj= 0
## Call: abund_hell ~ 1
##
## R2.adjusted
## <All variables> 0.039094
## + SmokingStatus 0.030093
## + Sex 0.016765
## + PackYears 0.013326
## <none> 0.000000
## + Age -0.002834
##
## Df AIC F Pr(>F)
## + SmokingStatus 1 -46.1 2.83 0.004 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Step: R2.adj= 0.03009
## Call: abund_hell ~ SmokingStatus
##
## R2.adjusted
## + Sex 0.04449
## <All variables> 0.03909
## <none> 0.03009
## + PackYears 0.02786
## + Age 0.02574La variable estado de fumador se selecciona mediante la selección de
avance utilizando la función ordiR2step(). Después de
ejecutar el procedimiento de selección hacia adelante, ajustamos el
modelo final parsimonioso de RDA y probamos con el modelo global, cada
eje y cada variable:
rda_final<- rda(abund_hell ~ SmokingStatus + Sex, data=throat_meta)anova(rda_final, step=1000)%>%
kbl(booktabs = TRUE, align = "c",
caption = "ANOVA test global")%>%
kable_styling(position = "center",
latex_options = c("hold_position", "striped"),
font_size = 12)| Df | Variance | F | Pr(>F) | |
|---|---|---|---|---|
| Model | 2 | 0.0356 | 2.374 | 0.001 |
| Residual | 57 | 0.4271 | NA | NA |
anova(rda_final, by="axis", step=1000)%>%
kbl(booktabs = TRUE, align = "c",
caption = "ANOVA test para ejes")%>%
kable_styling(position = "center",
latex_options = c("hold_position", "striped"),
font_size = 12)| Df | Variance | F | Pr(>F) | |
|---|---|---|---|---|
| RDA1 | 1 | 0.0226 | 3.010 | 0.002 |
| RDA2 | 1 | 0.0130 | 1.737 | 0.041 |
| Residual | 57 | 0.4271 | NA | NA |
anova(rda_final, by="terms", step=1000)%>%
kbl(booktabs = TRUE, align = "c",
caption = "ANOVA test para variables explicativas")%>%
kable_styling(position = "center",
latex_options = c("hold_position", "striped"),
font_size = 12)| Df | Variance | F | Pr(>F) | |
|---|---|---|---|---|
| SmokingStatus | 1 | 0.0215 | 2.873 | 0.001 |
| Sex | 1 | 0.0140 | 1.874 | 0.024 |
| Residual | 57 | 0.4271 | NA | NA |
Análisis de correspondencia restringido (CCA)
El CCA también se conoce como análisis de correspondencia canónica.
Desde su introducción en 1986, ha sido uno de los métodos de ordenación
más populares en ecología de comunidades y aceptado por los
investigadores del microbioma (ter Braak 1986). Al igual que el
RDA se relaciona con el PCA, el CAP se relaciona con el PCoA, y
el CCA se relaciona con el CA. CCA comparte las propiedades
básicas de CA y las combina en una ordenación restringida. CCA se
realiza mediante la función cca(). Su algoritmo se basa en
el de Legendre y Legendre (1998): preserva la distancia \(\chi^{2}\) entre muestras, y los taxones se
representan como puntos en los tripplots (Borcard et al. 2011). La
distancia \(\chi^{2}\) calculada se
somete a una regresión lineal ponderada sobre las variables de
restricción, y los valores ajustados se pasan al análisis de
correspondencia realizado mediante la descomposición de valores
singulares (svd). Por lo tanto, es una forma ponderada de RDA.
Al igual que RDA, hay dos tipos de sintaxis de CCA. Una es la sintaxis matricial simple (por defecto)
\[ cca(X,Y,Z) \]
donde \(X\) = matriz de datos de la comunidad o marco de datos, se debe dar; \(Y\) = matriz o marco de datos de restricción, normalmente de variables ambientales (puede omitirse); si la matriz \(Y\) se suministra, se utiliza para llevar a cabo el CCA, si no se suministra, se calcula el AC. \(Z\) = matriz condicionante o marco de datos (también puede omitirse). Si se suministra la matriz \(Z\), el efecto se partirá de la matriz de la comunidad de la comunidad.
Otra es la sintaxis de la fórmula:
\[ cca(formula, data, na.action = na.fail, subset = NULL) \]
donde, la fórmula es la fórmula típica del modelo: el lado izquierdo
de la fórmula debe ser la matriz de datos de la comunidad (\(X\)); el lado derecho define el modelo de
restricción: las variables de restricción pueden contener factores
ordenados o desordenados, interacciones entre las variables y funciones
de las variables; y las variables condicionantes pueden estar dadas
dentro de una condición de función especial para las variables
condicionantes (covariables) parcializadas antes del análisis. Así, los
siguientes comandos son equivalentes: \(cca(X,
Y, Z)\), \(cca(X ~ Y +
condition(Z))\), donde \(Y\) y
\(Z\) se refieren a las restricciones y
las condiciones respectivamente. Los datos son un marco de datos que
contiene las variables del lado derecho de la fórmula del modelo. La
función na.action() se utiliza para tratar los valores
perdidos en las restricciones o condiciones. El valor por defecto
(na.fail) es dejar de tener valores perdidos; na.omit es
eliminar todas las filas con valores perdidos; na.exclude
es mantener todas las observaciones pero dar NA para los
resultados que no se pueden calcular. Sin embargo, los valores perdidos
nunca se permiten en los datos de las comunidades dependientes. El
subconjunto se utiliza para subconjuntos de filas de datos. En esta
sección, utilizaremos los datos de los fumadores para realizar el CCA
mediante la función cca() del paquete vegan.
Hemos mencionado anteriormente que para realizar el CCA, la matriz de
variables ambientales debe ser suministrada, y de lo contrario la
función cca() calcula CA. Como recordamos, los datos de los
fumadores tienen dos conjuntos de datos: throat.otu.tab
(marco de datos de abundancia de la comunidad) y
throat.meta (metadatos que incluyen dos variables binarias
SmokingStatus y Sex, y dos variables continuas
Age y PackYears). Para ejecutar un CCA,
cargamos el paquete vegan y no transformamos los datos de
abundancia de la comunidad por el método Hellinger como hicimos en RDA.
De lo contrario, la distancia \(\chi^{2}\) no puede ser calculada y los
resultados no pueden ser interpretados.
En el siguiente CCA, los datos de abundancia de fumadores están
restringidos por cuatro variables ambientales ambientales en los datos
de throat.meta, incluyendo SmokingStatus,
Age, Sex, y PackYears.
smoker_cca <- cca(throat.otu.tab ~ ., throat_meta)
smoker_cca## Call: cca(formula = throat.otu.tab ~ SmokingStatus + Age + Sex +
## PackYears, data = throat_meta)
##
## Inertia Proportion Rank
## Total 4.9330 1.0000
## Constrained 0.3782 0.0767 4
## Unconstrained 4.5548 0.9233 55
## Inertia is scaled Chi-square
##
## Eigenvalues for constrained axes:
## CCA1 CCA2 CCA3 CCA4
## 0.1517 0.0912 0.0781 0.0572
##
## Eigenvalues for unconstrained axes:
## CA1 CA2 CA3 CA4 CA5 CA6 CA7 CA8
## 0.457 0.358 0.325 0.299 0.237 0.187 0.165 0.160
## (Showing 8 of 55 unconstrained eigenvalues)La primera sección después de la “Call:” da los coeficientes de contingencia media al cuadrado del análisis. En el CCA desempeñan el mismo papel que la varianza total en la RDA. Como la matriz de comunidad se convierte en distancia \(\chi^{2}\), las entradas son coeficientes de contingencia. La variación total antes de someter la matriz a regresión ponderada es de 4.933; ésta es la variación que podría explicarse. La variación en la matriz de la comunidad que se explica después de la regresión ponderada es 0.3782; esta es la variación que será explicada por los ejes en el CCA. La varianza de los residuos de la regresión es de 4.5548; se trata de la variación que no será explicada por los ejes en el CCA, que puede ser sometida a CA.
Aquí, observamos que el CCA no tiene mucho éxito porque sólo 0.0767 o
\(0.0767\) (véase la columna Proporción
y la fila Restringida) de la variación total de datos fue capturada en
el CCA por las cuatro variables. La varianza explicada por ejes
particulares puede encontrarse mediante la función
summary():
head(summary(smoker_cca))##
## Call:
## cca(formula = throat.otu.tab ~ SmokingStatus + Age + Sex + PackYears, data = throat_meta)
##
## Partitioning of scaled Chi-square:
## Inertia Proportion
## Total 4.933 1.0000
## Constrained 0.378 0.0767
## Unconstrained 4.555 0.9233
##
## Eigenvalues, and their contribution to the scaled Chi-square
##
## Importance of components:
## CCA1 CCA2 CCA3 CCA4 CA1 CA2 CA3 CA4
## Eigenvalue 0.1517 0.0912 0.0781 0.0572 0.4567 0.3581 0.3247 0.2989
## Proportion Explained 0.0308 0.0185 0.0158 0.0116 0.0926 0.0726 0.0658 0.0606
## Cumulative Proportion 0.0308 0.0492 0.0651 0.0767 0.1692 0.2418 0.3076 0.3682
## CA5 CA6 CA7 CA8 CA9 CA10 CA11 CA12
## Eigenvalue 0.2373 0.1865 0.1653 0.1604 0.1280 0.1207 0.1189 0.1101
## Proportion Explained 0.0481 0.0378 0.0335 0.0325 0.0259 0.0245 0.0241 0.0223
## Cumulative Proportion 0.4163 0.4542 0.4877 0.5202 0.5461 0.5706 0.5947 0.6170
## CA13 CA14 CA15 CA16 CA17 CA18 CA19 CA20
## Eigenvalue 0.1097 0.0999 0.0904 0.0859 0.0831 0.0783 0.0726 0.0665
## Proportion Explained 0.0222 0.0203 0.0183 0.0174 0.0168 0.0159 0.0147 0.0135
## Cumulative Proportion 0.6392 0.6595 0.6778 0.6952 0.7120 0.7279 0.7426 0.7561
## CA21 CA22 CA23 CA24 CA25 CA26 CA27
## Eigenvalue 0.0637 0.0612 0.0591 0.0555 0.0531 0.0499 0.04870
## Proportion Explained 0.0129 0.0124 0.0120 0.0113 0.0108 0.0101 0.00987
## Cumulative Proportion 0.7690 0.7814 0.7934 0.8047 0.8154 0.8256 0.83543
## CA28 CA29 CA30 CA31 CA32 CA33 CA34
## Eigenvalue 0.04757 0.04608 0.04326 0.04260 0.04030 0.03908 0.03808
## Proportion Explained 0.00964 0.00934 0.00877 0.00864 0.00817 0.00792 0.00772
## Cumulative Proportion 0.84507 0.85441 0.86318 0.87182 0.87999 0.88791 0.89563
## CA35 CA36 CA37 CA38 CA39 CA40 CA41
## Eigenvalue 0.03580 0.03498 0.0340 0.0331 0.03105 0.03089 0.02951
## Proportion Explained 0.00726 0.00709 0.0069 0.0067 0.00629 0.00626 0.00598
## Cumulative Proportion 0.90289 0.90998 0.9169 0.9236 0.92987 0.93613 0.94211
## CA42 CA43 CA44 CA45 CA46 CA47 CA48
## Eigenvalue 0.02793 0.02646 0.02574 0.02359 0.02256 0.02217 0.02092
## Proportion Explained 0.00566 0.00536 0.00522 0.00478 0.00457 0.00449 0.00424
## Cumulative Proportion 0.94778 0.95314 0.95836 0.96314 0.96771 0.97221 0.97645
## CA49 CA50 CA51 CA52 CA53 CA54 CA55
## Eigenvalue 0.02013 0.0192 0.01798 0.01757 0.01556 0.01342 0.01228
## Proportion Explained 0.00408 0.0039 0.00365 0.00356 0.00316 0.00272 0.00249
## Cumulative Proportion 0.98053 0.9844 0.98807 0.99164 0.99479 0.99751 1.00000
##
## Accumulated constrained eigenvalues
## Importance of components:
## CCA1 CCA2 CCA3 CCA4
## Eigenvalue 0.152 0.0912 0.0781 0.0572
## Proportion Explained 0.401 0.2411 0.2065 0.1512
## Cumulative Proportion 0.401 0.6422 0.8488 1.0000
##
## Scaling 2 for species and site scores
## * Species are scaled proportional to eigenvalues
## * Sites are unscaled: weighted dispersion equal on all dimensions
##
##
## Species scores
##
## CCA1 CCA2 CCA3 CCA4 CA1 CA2
## 4695 0.0116 1.177 0.475 0.1308 0.0491 0.0182
## 2983 0.5220 0.510 1.025 1.2222 0.1317 0.6923
## 2554 0.0277 1.202 0.666 0.0181 -0.1559 1.0291
## 3315 -0.0679 0.382 0.474 -0.1129 -0.4215 3.1688
## 879 -1.1588 0.745 -0.160 -0.8216 -0.6967 -0.5789
## 1313 -0.0124 -1.048 -0.370 0.0333 -0.2758 3.2476
## ....
##
##
## Site scores (weighted averages of species scores)
##
## CCA1 CCA2 CCA3 CCA4 CA1 CA2
## ESC_1.1_OPL -2.510 2.1425 0.415 -2.318203 2.6670 -0.6690
## ESC_1.3_OPL 1.309 -0.0894 -0.376 -0.176963 0.0602 -0.3029
## ESC_1.4_OPL 1.445 -0.0145 -0.364 0.078627 0.2178 -0.6817
## ESC_1.5_OPL 1.538 0.4089 -0.589 0.896660 -0.4024 -1.2717
## ESC_1.6_OPL 2.579 0.8192 -2.732 0.254050 0.2062 -0.5353
## ESC_1.10_OPL 0.385 -0.9651 -0.891 0.000538 -0.1476 -0.0726
## ....
##
##
## Site constraints (linear combinations of constraining variables)
##
## CCA1 CCA2 CCA3 CCA4 CA1 CA2
## ESC_1.1_OPL -1.125 0.7768 -0.132 -0.961 2.6670 -0.6690
## ESC_1.3_OPL 1.123 0.0630 0.341 0.414 0.0602 -0.3029
## ESC_1.4_OPL 1.067 0.0273 0.578 0.567 0.2178 -0.6817
## ESC_1.5_OPL -0.029 1.8898 0.889 0.728 -0.4024 -1.2717
## ESC_1.6_OPL 1.261 0.1173 -0.834 0.195 0.2062 -0.5353
## ESC_1.10_OPL 1.598 0.2708 -3.324 -0.438 -0.1476 -0.0726
## ....
##
##
## Biplot scores for constraining variables
##
## CCA1 CCA2 CCA3 CCA4 CA1 CA2
## SmokingStatusSmoker 0.831 0.452 0.0592 0.3183 0 0
## Age 0.202 -0.178 -0.7739 0.5730 0 0
## SexMale 0.682 -0.730 0.0444 0.0112 0 0
## PackYears 0.615 0.163 -0.7705 0.0307 0 0
##
##
## Centroids for factor constraints
##
## CCA1 CCA2 CCA3 CCA4 CA1 CA2
## SmokingStatusNonSmoker -0.746 -0.406 -0.0531 -0.28569 0 0
## SmokingStatusSmoker 0.926 0.503 0.0660 0.35459 0 0
## SexFemale -0.922 0.985 -0.0599 -0.01516 0 0
## SexMale 0.505 -0.540 0.0328 0.00831 0 0La varianza explicada por los cuatro primeros ejes restringidos (CCA1, CCA2, CCA3, CCA4): 3,07, 1,85 1,58 y 1,16%, respectivamente. La varianza explicada por el primer eje no restringido (CCA1) es del 9,26%. La sección de “Valores propios restringidos acumulados” proporciona los valores propios asociados a la proyección. Como tenemos cuatro variables en nuestro marco de datos ambiental, hay cuatro valores propios limitados. Como se muestra aquí, el primer eje representa aproximadamente el 40% de la variación restringida, con el segundo eje el segundo eje con un 24%, el tercero con un 21% y el cuarto con un 15%.
A continuación, trazamos los resultados del CCA para el escalado de
tipo 1 utilizando la función plot(). En el gráfico,
queremos mostrar las restricciones lineales o “puntuaciones LC”
(display = "lc") y los centroides de los niveles de las
variables factoriales (display = "cn").
plot(smoker_cca, scaling=1, display = c("lc","cn"), main="Biplot CCA - scaling 1") 
CCA biplot de la abundancia de la base smoker.throat restringido por las variables SmokingStatus, Age, Sex, PackYear con escalado tipo 1
La escala de tipo 1 muestra cuatro grupos de muestras, con el grupo de fumadores vinculado a PackYears. En este análisis, el primer eje se asocia con el aumento de PackYears mientras que el segundo está asociado con la disminución de la Edad.
Al igual que en RDA, realizamos una prueba de permutación para el modelo global, cada eje y cada variable explicativa. Para asegurarse de obtener el mismo resultado de la prueba de permutación cada vez que se ejecute, establezca la misma semilla.
set.seed (123)
anova(smoker_cca, step=1000) %>%
kbl(booktabs = TRUE, align = "c",
caption = "ANOVA test global: El resultado no es estadísticamente significativo")%>%
kable_styling(position = "center",
latex_options = c("hold_position", "striped"),
font_size = 12)| Df | ChiSquare | F | Pr(>F) | |
|---|---|---|---|---|
| Model | 4 | 0.3782 | 1.142 | 0.14 |
| Residual | 55 | 4.5548 | NA | NA |
Como el resultado anterior no es estadísticamente significativo, entonces realizamos una prueba de significancia por cada eje.
set.seed (123)
anova (smoker_cca, by = 'axis',step=1000) %>%
kbl(booktabs = TRUE, align = "c",
caption = "ANOVA test por ejes: El primer eje es significativo")%>%
kable_styling(position = "center",
latex_options = c("hold_position", "striped"),
font_size = 12)| Df | ChiSquare | F | Pr(>F) | |
|---|---|---|---|---|
| CCA1 | 1 | 0.1517 | 1.8318 | 0.093 |
| CCA2 | 1 | 0.0912 | 1.1009 | 0.824 |
| CCA3 | 1 | 0.0781 | 0.9432 | 0.853 |
| CCA4 | 1 | 0.0572 | 0.6906 | 0.941 |
| Residual | 55 | 4.5548 | NA | NA |
Ahora probamos la signifiancia de cada variable explicativa.
set.seed(123)
anova (smoker_cca, by = 'terms')%>%
kbl(booktabs = TRUE, align = "c",
caption = "ANOVA test por variables explicativass: SmokingStatus
es significativo con un p-value de 0.009")%>%
kable_styling(position = "center",
latex_options = c("hold_position", "striped"),
font_size = 12)| Df | ChiSquare | F | Pr(>F) | |
|---|---|---|---|---|
| SmokingStatus | 1 | 0.1295 | 1.5639 | 0.008 |
| Age | 1 | 0.0737 | 0.8895 | 0.676 |
| Sex | 1 | 0.1049 | 1.2663 | 0.078 |
| PackYears | 1 | 0.0701 | 0.8468 | 0.700 |
| Residual | 55 | 4.5548 | NA | NA |
Por último, vamos a hacer una selección hacia delante utilizando la
función ordistep() del paquete vegan. Como
dijimos en RDA, la función ordistep() es aplicable con las
funciones rda(), cca() o
cmdscale(). La comparación de variables se basa en el
criterio AIC y en los valores p de la prueba de permutación de Monte
Carlo.
ordistep(cca(throat.otu.tab ~ 1, data=throat_meta), scope=formula(smoker_cca),
direction="forward", pstep=1000)##
## Start: throat.otu.tab ~ 1
##
## Df AIC F Pr(>F)
## + SmokingStatus 1 539 1.56 0.025 *
## + Sex 1 539 1.44 0.030 *
## + PackYears 1 539 1.28 0.190
## + Age 1 540 0.89 0.585
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Step: throat.otu.tab ~ SmokingStatus
##
## Df AIC F Pr(>F)
## + Sex 1 540 1.28 0.06 .
## + PackYears 1 540 0.97 0.51
## + Age 1 540 0.89 0.66
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1## Call: cca(formula = throat.otu.tab ~ SmokingStatus, data =
## throat_meta)
##
## Inertia Proportion Rank
## Total 4.9330 1.0000
## Constrained 0.1295 0.0263 1
## Unconstrained 4.8035 0.9737 58
## Inertia is scaled Chi-square
##
## Eigenvalues for constrained axes:
## CCA1
## 0.1295
##
## Eigenvalues for unconstrained axes:
## CA1 CA2 CA3 CA4 CA5 CA6 CA7 CA8
## 0.484 0.383 0.332 0.302 0.249 0.188 0.177 0.161
## (Showing 8 of 58 unconstrained eigenvalues)smoker_cca_final <- cca(throat.otu.tab ~ SmokingStatus, data=throat_meta)
anova.cca(smoker_cca_final, step=1000)%>%
kbl(booktabs = TRUE, align = "c",
caption = "ANOVA test global")%>%
kable_styling(position = "center",
latex_options = c("hold_position", "striped"),
font_size = 12)| Df | ChiSquare | F | Pr(>F) | |
|---|---|---|---|---|
| Model | 1 | 0.1295 | 1.564 | 0.004 |
| Residual | 58 | 4.8035 | NA | NA |
anova.cca(smoker_cca_final, step=1000, by="terms") %>%
kbl(booktabs = TRUE, align = "c",
caption = "ANOVA test por variables explicativass")%>%
kable_styling(position = "center",
latex_options = c("hold_position", "striped"),
font_size = 12)| Df | ChiSquare | F | Pr(>F) | |
|---|---|---|---|---|
| SmokingStatus | 1 | 0.1295 | 1.564 | 0.004 |
| Residual | 58 | 4.8035 | NA | NA |
La condición de fumador alcanza el mismo nivel de significación con un \(p-\)valorde \(0,007\). Sin embargo, el modelo de parsimonia ha dado sus frutos con un residuo mayor.
Análisis restringido de coordenadas principales (CAP)
CAP (también llamado análisis restringido de proximidades en el
paquete vegan), es un método de ordenación similar a RDA.
Es simplemente un análisis de redundancia de los resultados del
análisis de coordenadas principales (o escalamiento multidimensional
métrico) (Anderson y Willis 2003). CAP permite índices de
disimilitud no euclidianos, como la distancia Manhattan o Bray-Curtis.
Si se especifica la distancia euclidiana como método de ordenación, los
resultados serán idénticos a los de RDA.
La función capscale() del paquete vegan se
utiliza para implementar CAP. Necesita una matriz de disimilitud como
conjunto de datos de entrada, que puede calcularse mediante las
funciones vegdist(), dist(), o cualquier otro
método que produzca matrices de similitud. Se necesitan dos pasos
necesarios: primero, utilizar la función cmdscale() para
ordenar la matriz de disimilitud, luego utilizar RDA para analizar estos
resultados. A diferencia de RDA, en el que se pueden utilizar tanto la
matriz como la la función capscale() sólo se puede llamar
con la sintaxis de la fórmula. A continuación se muestra un uso de la
función capscale().
capscale(formula, data, distance = "bray", dfun = vegdist)
donde, formula es una fórmula típica del modelo como se define en
rda() y cca(). El lado izquierdo de la fórmula
debe ser una matriz de datos de la comunidad (marco) o una matriz de
disimilitud que puede estimarse con la función vegdist() o
dist(). Si el lado izquierdo de la fórmula es una matriz de
datos (datframe) en lugar de una matriz de disimilitud, entonces se debe
proporcionar un índice de disimilitud (o distancia). El lado derecho de
la fórmula define las restricciones. Las variables de restricción pueden
ser variables continuas, factores, términos de interacción o una
condición de término especial utilizada para definir las variables que
se van a particionar. Los datos son un data frame que contiene las
variables del lado derecho de la fórmula del modelo. La función
dfun es la función de distancia o función de disimilitud
utilizada.
El CAP básico se puede realizar con los siguientes códigos. Las
variables de restricción incluyen dos variables binarias
(SmokingStatus y Sex), y dos variables
continuas (PackYears y Age). La condición de
término especial (Age) se utiliza para particionar el
efecto de la edad.
throat_meta <- select(throat.meta, SmokingStatus, Age, Sex, PackYears)
throat_cap <- capscale(throat.otu.tab ~ SmokingStatus + Sex + PackYears + Condition(Age), throat_meta,
dist="bray")
throat_cap## Call: capscale(formula = throat.otu.tab ~ SmokingStatus + Sex +
## PackYears + Condition(Age), data = throat_meta, distance = "bray")
##
## Inertia Proportion Rank
## Total 14.0932 1.0000
## Conditional 0.2084 0.0148 1
## Constrained 1.2623 0.0896 3
## Unconstrained 13.0648 0.9270 46
## Imaginary -0.4425 -0.0314 13
## Inertia is squared Bray distance
## Species scores projected from 'throat.otu.tab'
##
## Eigenvalues for constrained axes:
## CAP1 CAP2 CAP3
## 0.731 0.376 0.155
##
## Eigenvalues for unconstrained axes:
## MDS1 MDS2 MDS3 MDS4 MDS5 MDS6 MDS7 MDS8
## 2.278 1.994 1.144 0.925 0.749 0.615 0.534 0.473
## (Showing 8 of 46 unconstrained eigenvalues)Los tres primeros ejes se denominan “CAP1”, “CAP2” y “CAP3”, y el
siguiente es el MDS original. Podemos ver que las tres variables
restringidas (SmokingStatus, Sex y
PackYears) explican el 8,83% de la variación total del
conjunto de datos.
set.seed(123)
anova1 <- anova(throat_cap)
p_valor1 <- anova1$`Pr(>F)`
anova1 %>% kbl(booktabs = TRUE, align = "c",
caption = "ANOVA test")%>%
kable_styling(position = "center",
latex_options = c("hold_position", "striped"),
font_size = 12)| Df | SumOfSqs | F | Pr(>F) | |
|---|---|---|---|---|
| Model | 3 | 1.262 | 1.771 | 0.007 |
| Residual | 55 | 13.065 | NA | NA |
El modelo es estadísticamente significativo con un \(p-\)valor de 0.007, NA. La función
plot(), genera la siguiente figura, pero no es
informativa
plot(throat_cap)
Plot básico de CAP en los datos de fumadores.
Como estamos interesados en las diferentes disimilitudes entre los fumadores y los no fumadores, vamos a extraer la información del grupo.
groups <- throat.meta$SmokingStatus
groups## [1] NonSmoker Smoker Smoker Smoker Smoker Smoker NonSmoker
## [8] NonSmoker NonSmoker NonSmoker Smoker NonSmoker NonSmoker Smoker
## [15] NonSmoker Smoker Smoker NonSmoker NonSmoker NonSmoker NonSmoker
## [22] NonSmoker NonSmoker NonSmoker NonSmoker NonSmoker NonSmoker NonSmoker
## [29] NonSmoker NonSmoker NonSmoker NonSmoker NonSmoker Smoker NonSmoker
## [36] NonSmoker NonSmoker NonSmoker NonSmoker Smoker NonSmoker Smoker
## [43] NonSmoker Smoker Smoker Smoker Smoker Smoker Smoker
## [50] Smoker Smoker Smoker Smoker Smoker Smoker Smoker
## [57] Smoker NonSmoker Smoker Smoker
## Levels: NonSmoker SmokerEn el siguiente grupo de códigos R, la función plot()
genera un diagrama de ordenación CAP vacío; la función
points() añade puntos al diagrama de ordenación (función de
trazado de bajo nivel) creado por la función plot. La función
ordispider() crea spiderplot conectando los miembros
individuales del grupo con el centroide del grupo. La función
ordiellipse() rodea las nubes de puntos dentro del grupo
mediante una elipse como las envolventes.
# Correr todo junto
plot(throat_cap, type="n")
points(throat_cap, col=as.numeric(as.factor(groups)),
pch=as.numeric(as.factor(groups)))
ordispider(throat_cap, groups, lty=2, col="grey", label=T)
ordiellipse(throat_cap, groups, lty=2, col="grey", label=F) 
Al igual que en RDA y CCA, vamos a realizar una prueba de permutación para el modelo global global, cada eje y cada variable explicativa. Para asegurarse de obtener el mismo resultado de la prueba de permutación cada vez que se ejecute, establezca la misma semilla.
set.seed (123)
anova_global <- anova(throat_cap, step=1000)
pvalue_global <- anova_global$`Pr(>F)`
anova_global %>% kbl(booktabs = TRUE, align = "c",
caption = "ANOVA test global")%>%
kable_styling(position = "center",
latex_options = c("hold_position", "striped"),
font_size = 12)| Df | SumOfSqs | F | Pr(>F) | |
|---|---|---|---|---|
| Model | 3 | 1.262 | 1.771 | 0.007 |
| Residual | 55 | 13.065 | NA | NA |
Observamos que el resultado es estadísticamente significativocon un \(p-\)valor de 0.007, NA. Entonces probamos significancia con cada eje.
set.seed (123)
anova_ejes <- anova(throat_cap, by="axis", step=1000)
pvalue_CA1 <- anova_ejes$`Pr(>F)`[1]
pvalue_CA2 <- anova_ejes$`Pr(>F)`[2]
anova_ejes %>%
kbl(booktabs = TRUE, align = "c",
caption = "ANOVA test por ejes")%>%
kable_styling(position = "center",
latex_options = c("hold_position", "striped"),
font_size = 12)| Df | SumOfSqs | F | Pr(>F) | |
|---|---|---|---|---|
| CAP1 | 1 | 0.7308 | 3.0763 | 0.014 |
| CAP2 | 1 | 0.3764 | 1.5844 | 0.236 |
| CAP3 | 1 | 0.1552 | 0.6535 | 0.871 |
| Residual | 55 | 13.0648 | NA | NA |
Hay significancia en el primer eje con un \(p-\)valor de 0.014 y es margialmente significativo en el segundo eje con un \(p-\)valor de 0.236. Ahora probamos con cada valiable explicativa.
set.seed (123)
anova_varibales <- anova(throat_cap, by="terms", step=1000)
pvalue_SmockingStatus <- anova_varibales$`Pr(>F)`[1]
pvalue_Sex <- anova_varibales$`Pr(>F)`[2]
anova_varibales %>%
kbl(booktabs = TRUE, align = "c",
caption = "ANOVA test por variables explicativas")%>%
kable_styling(position = "center",
latex_options = c("hold_position", "striped"),
font_size = 12) | Df | SumOfSqs | F | Pr(>F) | |
|---|---|---|---|---|
| SmokingStatus | 1 | 0.6305 | 2.6544 | 0.002 |
| Sex | 1 | 0.4438 | 1.8683 | 0.033 |
| PackYears | 1 | 0.1880 | 0.7914 | 0.697 |
| Residual | 55 | 13.0648 | NA | NA |
Podemos ver que tanto SmokingStatus como
Sex son significativos con \(p-\)valores de 0.002 y 0.033
respectivamente.
Ahora, hagamos la selección hacia delante utilizando la función
ordistep() del paquete vegan. Como dijimos en
RDA y CCA, la función compara las variables basándose en el criterio AIC
y los \(p-\)valores de la prueba de
permutación de Monte Carlo.
step_forward <- ordistep(capscale(throat.otu.tab ~ 1, data=throat_meta),
scope=formula(throat_cap), direction="forward", pstep=1000)##
## Start: throat.otu.tab ~ 1
##
## Df AIC F Pr(>F)
## + SmokingStatus 1 712 2.05 0.06 .
## + Sex 1 712 1.68 0.09 .
## + PackYears 1 712 1.27 0.26
## + Condition(Age) 1 713 0.00
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# Seleccionamos la variable de Smoking status, el modelo final es ajustado como sigue
cap_final<- capscale(throat.otu.tab ~ SmokingStatus, throat_meta, dist="bray")
set.seed(123)
anova_M <- anova(cap_final, step=1000)
p_valueM <- anova_M$`Pr(>F)`
anova_M %>%
kbl(booktabs = TRUE, align = "c",
caption = "ANOVA test")%>%
kable_styling(position = "center",
latex_options = c("hold_position", "striped"),
font_size = 12) | Df | SumOfSqs | F | Pr(>F) | |
|---|---|---|---|---|
| Model | 1 | 0.6533 | 2.729 | 0.002 |
| Residual | 58 | 13.8823 | NA | NA |
El estatus de fumador alcanza la significación con un \(p-\)valor de 0.002, NA, pero el modelo de parsimonia ha dado sus frutos con un residuo mayor.
Resumen y discusión
Utilizamos conjuntos de datos de ratones y humanos para ilustrar el
análisis exploratorio de los datos del microbioma. En primer lugar,
utilizamos el paquete phyloseq para ilustrar los cinco
gráficos, incluyendo la riqueza, la barra de abundancia, el mapa de
calor, la red y el árbol filogenético. A continuación, introdujimos
varias familias de métodos de agrupación disponibles en los
estudios de ecología y microbioma, y nos centramos en ilustrar cuatro
métodos de agrupación (aglomerativo de enlace único,
aglomerativo de enlace completo, aglomerativo de tinta
media y aglomerativo de varianza mínima de Ward). También
describimos brevemente la relación entre clustering, ordenación y medida
de distancia. Por último, ilustramos las ordenaciones más
comunes sin restricciones y con restricciones:
PCA, PCoA, NMDS,
CA, RDA, CCA y
CAP. Las características de las ordenaciones no
restringidas se consideran “exploratorias”. Sin embargo, las capacidades
de las ordenaciones restringidas van más allá del mero análisis
exploratorio de datos, y se convierten también en pruebas de hipótesis.
Hemos introducido las diferentes características de las ordenaciones no
restringidas y restringidas.
Ilustramos el análisis exploratorio y la prueba de hipótesis (en
términos de ordenaciones restringidas) de los datos del microbioma
mediante el uso del paquete para analizar los datos del censo del
microbioma phyloseq y los paquetes que comúnmente se
utilizan en ecología, como el paquete vegan. Las funciones
de extensión, es decir cleanplot.pca() y
evplot() son atractivas y dos criterios para evaluar el
rendimiento de las ordenaciones (el criterio de Kaiser-Guttman
y el modelo broken-stick) son útiles para evaluar los análisis
de ordenación.